Presentation Title

For a bounded Lipschitz domain $\Omega$ in $\R^d$ with $d \in \lbrace 2,3 \rbrace$, we consider \begin{align} \partial_t u(t,x) -\Delta_x u(t,x) &= f(t,x), \quad (t,x) \in Q = I \times \Omega, \nonumber \\ u(t,x) &= 0, \quad \quad \quad \, (t,x) \in \Sigma = (0,T) \times \partial\Omega, \nonumber \\ u(0,x) &= 0, \quad \quad \quad \, \, x \in \Omega. \end{align} Discretization in space first leads to ODE system \begin{align} \mathbf{M} \frac{d\mathbf{u}}{dt}(t) + \mathbf{K} \mathbf{u}(t) = \mathbf{f}(t), \quad \mathbf{u}(0) = 0, \end{align} where $\mathbf{M},\mathbf{K}$ are matrices.

• Classical ODE solvers are inherently sequential.
• We aim to find a space-time discretization

For a bounded Lipschitz domain $\Omega$ in $\R^d$ with $d \in \lbrace 2,3 \rbrace$, we consider \begin{align} \partial_t u(t,x) -\Delta_x u(t,x) &= f(t,x), \quad (t,x) \in Q = I \times \Omega, \nonumber \\ u(t,x) &= 0, \quad \quad \quad \, (t,x) \in \Sigma = (0,T) \times \partial\Omega, \nonumber \\ u(0,x) &= 0, \quad \quad \quad \, \, x \in \Omega. \end{align} Eigenfunction decomposition $u(t,x)= \sum_{n=0}^{\infty} u_n(t) \varphi_n(x)$ with \begin{align} -\Delta_x \varphi_n = \mu_n \varphi_n, \quad \varphi \mid_{\partial \Omega} = 0 \end{align} yields mode-wise equation \begin{align} \partial_t u_n(t) + \mu_n u_n(t) = f_n(t), \quad u_n(0)=0. \end{align}

• Our approach is not limited to the heat equation.

• \(\mu>0 \) represents spatial operator.
• Exact solution \(u(t) = \int_0^t e^{-\mu (t-s)} f(s) \, ds \).
• Variational setting?

Primal Formulation: Dual Formulation:

• Is there something in between? \(\Longrightarrow\) Lions-Magenes setting^1,2

Use the Fourier series $ u(t) = \sum_{k=0}^\infty u_k \sin\left(\left(k \pi+ \frac{\pi}{2}\right)\frac{t}{T}\right)$ to define for \(s \geq 0\) the mapping \[ \begin{align} \mathcal{H}_T \, : \, H^s_{0,}(I) &\to H^s_{,0}(I), \\ u &\mapsto \sum_{k=0}^\infty u_k \cos\left(\left(k \pi+ \frac{\pi}{2}\right)\frac{t}{T}\right) \end{align} \]

• \(\mathcal{H}_T\) is an isomorphism.
• $\lVert u \rVert_{H^{1/2}_{0,}(I)}^2 = \langle \partial_t u, \mathcal{H}_T u \rangle_I = \langle \mathcal{H}_T u, \partial_t u \rangle_I.$

Find \(u \in X = H^{1/2}_{0,}(I)\) such that

• Elliptic: \( b(u,u) \geq \lVert u \rVert_X^2 \)
• Bounded: \( b(u,v) \leq C (1 + \mu) \lVert u \rVert_X \lVert v \rVert_X \)
• Boundedness depends on $\mu$, but ellipticity does not!
  \(\Longrightarrow\) problem in the PDE case! (norm-gap)
  \(\Longrightarrow\) can we find a better transformation?

Consider the operator \[ \begin{align} T_w \, : \, \C & \to \C, \\ z & \mapsto wz \,. \end{align} \] Variational formulation to find \(u \in \C\) such that \[ b(u,v) = \langle T_w u, v \rangle_\C = \langle f,v \rangle_\C \quad \forall v \in \C. \] with $\langle v, w \rangle_\C = v \overline{w}$.

Note, that \( b(u,u) = w \lvert u \rvert^2 \) is indefinite.
\(\Longrightarrow\) No Lax-Milgram!

• Can we find a transformation such that ellipticity is recovered?

Using $w = e^{i \varphi} \lvert w \rvert$, we split $T_w$ into \[ \begin{align} T_w = U_w \lvert T_w \rvert \end{align} \] with \[ \begin{align} U_w v = e^{i \varphi} v, \quad \lvert T_w \rvert v = \lvert w \rvert \, v. \end{align} \]

• $U^*v = e^{-i \varphi} v$
• $U^* U = I \quad $ $\Longrightarrow$ $U$ is unitary
• $\lVert U v \rVert_\C^2 = \langle U v, U v \rangle_\C = \langle v, v \rangle_\C = \lVert v \rVert_\C^2 \quad$ $\Longrightarrow$ $U$ is an isometry

Since $U_w$ is unitary, we can augment $b(\cdot,\cdot)$ to \[ \begin{align} \tilde{b}(u,v) &= \langle T_w u, U_w v \rangle_\C = \langle U_w^* T_w u, v \rangle_\C= \langle U_w^* U_w \lvert T_w \rvert u, v \rangle_\C =\langle \lvert T_w \rvert u, v \rangle_\C \quad \forall u,v \in \C \end{align} \] and obtain ellipticity: \[ \tilde{b}(u,u) = \lVert w \rVert_\C \lVert u \rVert_\C^2. \] The variational problem to find \(u \in \C\) such that \[ \tilde{b}(u,v) = \langle f, U_w v \rangle_\C \quad \forall v \in \C \] is well-posed by Lax-Milgram. $\Longrightarrow$ Can we use a similar approach for our surrogate ODE?

Let $T \, : \, \mathcal{H} \supset \dom(T) \to \mathcal{H}$ be a closed and densely defined operator. Then, there exists
a positive self adjoint operator $\lvert T \rvert$, i.e. \begin{align} \langle \lvert T \rvert u, u \rangle_{\mathcal{H}} > 0 \quad \forall u \in \dom(\lvert T \rvert) \setminus \lbrace 0 \rbrace, \end{align} with domain $\dom(\lvert T \rvert) = \dom(T)$ and a partial isometry $U$ such that \begin{align} T = U \lvert T \rvert, \end{align} where the initial space of $U$ is $\left(\ker(T)\right)^\perp$ and the final space is $\overline{\ran}(T)$.
The absolute value operator is given by \begin{align} \lvert T \rvert = \left(T^* T\right)^{1/2}. \end{align}

The ODE operator acts on $L^2(I)$ as \begin{align} B_\mu \, : \, L^2(I) \supset \dom(B_\mu) &\to L^2(I), \\ u &\mapsto \partial_t u + \mu u, \end{align} with $\dom(B_\mu) = \left\lbrace v \in L^2(I) \, : \, \partial_t v \in L^2(I), \, v(0) = 0 \right\rbrace = H^1_{0,}(I)$. Its adjoint (in the sense of unbounded operators) has domain $\dom(B_\mu^*) = H^1_{,0}(I)$ and is given by \begin{align} B_\mu^* \, : \, L^2(I) \supset \dom(B_\mu^*) &\to L^2(I), \\ v &\mapsto -\partial_t v + \mu v. \end{align}

• $B_\mu$ and $B_\mu^*$ are closed and densely defined operators on $L^2(I)$.
• Both are isomorphisms from their domain to $L^2(I)$. $\Longrightarrow U$ is unitary.

\begin{align} S_\mu = B_\mu^* B_\mu = \left(-\partial_t + \mu \right)\left(\partial_t + \mu \right) = -\partial_{tt} + \mu^2, \end{align} with domain \begin{align} \dom(S_\mu) &= \left\lbrace v \in \dom(B_\mu) \, : \, B_\mu v \in \dom(B_\mu^*) \right\rbrace \\ &= \left \lbrace v \in H^1_{0,}(I) \, : \, \partial_t v + \mu v \in H^1_{,0}(I) \right\rbrace \\ &= \left\lbrace v \in H^2(I) \, : \, v(0) = v'(T) + \mu v(T) = 0 \right\rbrace. \end{align} Eigenvalue problem with Robin type boundary conditions: \begin{align} -\partial_{tt} \phi + \mu^2 \phi &= \lambda \phi, \\ \phi(0) = \phi'(T) + \mu \phi(T) &= 0. \end{align}

Solution of the eigenvalue problem: \begin{align} \phi_k(t) &= c_k \sin(\omega_k t), \quad \tan(\omega_k T) = -\frac{\omega_k}{\mu}, \\ \lambda_k &= \omega_k^2 + \mu^2, \qquad k = 0,1,2,\ldots \end{align}

\begin{align} \lvert B_\mu \rvert \sin(\omega_k t) &= \sqrt{{\omega_k^2 + \mu^2}} \, \sin(\omega_k t), \\ B_\mu \sin(\omega_k t) &= \omega_k \cos(\omega_k t) + \mu \sin(\omega_k t), \\ \Longrightarrow U_\mu \sin(\omega_k t) &= \frac{\omega_k}{\sqrt{\omega_k^2 + \mu^2}} \cos(\omega_k t) + \frac{\mu}{\sqrt{\omega_k^2 + \mu^2}} \sin(\omega_k t) \end{align} Variational formulation to find $u \in X=\dom(\lvert B_\mu \rvert^{1/2})$ such that \begin{align} b(u,v) = \langle B_\mu u, U_\mu v \rangle_I = \langle f, U_\mu v \rangle_I \quad \forall v \in X. \end{align}

• $b(u,v) = \langle \lvert B_\mu \rvert u, v \rangle_I = \langle \lvert B_\mu \rvert^{1/2} u, \lvert B_\mu \rvert^{1/2} v \rangle_I$ is symmetric, bounded and elliptic.
• $\lVert u \rVert_X^2 \simeq \lVert u \rVert_{H^{1/2}_{0,}(I)}^2 + \mu \lVert u \rVert_{L^2(I)}^2$ includes $\mu$.

Introduce a discrete variational problem to find $u \in X_s$ such that \begin{align} b_h(u,v) = \langle \partial_t u + \mu u, (I + \mathcal{H}_T) v \rangle_I = \langle f, (I + \mathcal{H}_T) v \rangle_I \quad \forall v \in X. \end{align}

• Well-defined for $X_s = H^{1/2+\varepsilon}_{0,}(I)$, $X=H^{1/2}_{0,}(I)$ with any $\varepsilon > 0$.
• Elliptic in full norm: $b(u,u) \geq \lVert u \rVert_X^2 = \lVert u \rVert_{H^{1/2}_{0,}(I)}^2 + \mu \lVert u \rVert_{L^2(I)}^2$.
• Galerkin discretization $u_h, v_h \in X_h = S_{h_t; 0,}^{p_t}(I) \cap X \subset X_s$ satisfies \begin{align} \lVert u - u_h \rVert_{X}\leq C \sqrt{1+h_t\mu} \, h_t^{\min\lbrace r, p_t + 1 \rbrace - 1/2} \lVert u \rVert_{H^r(I, \R)}, \end{align} and $b_h(u_h, v_h) \leq C \lVert u_h \rVert_X \lVert v_h \rVert_X$ by inverse inequality.

Use standard continuous piecewise polynomial spaces $X_h = S_{h_t; 0,}^{p_t}(I) \cap X = \span\lbrace \varphi_i \rbrace_{i=1}^{N_T}$ to obtain the equivalent linear system \begin{align} \mathbf{B} \mathbf{u} = \mathbf{f}, \end{align} \begin{align} (\mathbf{B})_{i,j} &= b_h(\varphi_j, \varphi_i) = \langle \partial_t \varphi_j + \mu \varphi_j, (I + \mathcal{H}_T) \varphi_i \rangle_I, \\ (\mathbf{f})_i &= \langle f, (I + \mathcal{H}_T) \varphi_i \rangle_I. \end{align}

• Matrices involving $\mathcal{H}_T$ are dense, but on regular grids FFT techniques reduce complexity to $\mathcal{O}(N_T \log(N_T))$. [pip: hmodfft]
• Fast approximation also works for non-linear problems and varying coefficients.
• Fast matrix-vector multiplication is possible, so preconditioning becomes the next step.

Boundedness and ellipticity in the $X$-norm imply \begin{align} \lVert u_h \rVert_{H^{1/2}_{0,}(I)}^2 + \mu\lVert u_h \rVert_{L^2(I)}^2 \leq \langle \mathbf{B} \mathbf{u}, \mathbf{u} \rangle_2 \leq \lVert u_h \rVert_{H^{1/2}_{0,}(I)}^2 + \mu \lVert u_h \rVert_{L^2(I)}^2. \end{align} Precondition with the inverse of the norm-inducing matrix $ \mathbf{D}_\mu = \mathbf{A}_t + \mu \mathbf{M}_t$ \begin{align} (\mathbf{A}_t)_{i,j} = \langle \partial_t \varphi_j, \mathcal{H}_T \varphi_i \rangle_I, \quad (\mathbf{M}_t)_{i,j} = \langle \varphi_j, \varphi_i \rangle_I. \end{align} We can prove that the symmetrically scaled system $ \mathbf{\tilde{B}} \, \mathbf{w} = \mathbf{g}$, \begin{align} \mathbf{\tilde{B}} = \mathbf{D}_\mu^{-1/2} \, \mathbf{B} \, \mathbf{D}_\mu^{-1/2}, \quad \mathbf{w} = \mathbf{D}_\mu^{1/2} \, \mathbf{u}, \quad \mathbf{g} = \mathbf{D}_\mu^{-1/2} \, \mathbf{f} \end{align} admits a bounded number of GMRES iterations.

Left preconditioning works well in practice. Behavior for $N_t = 200$ and $p_t = 1$ for different values of $\mu$ (right-hand side set to $1$):

Left preconditioning works well in practice. Replace $\mathbf{D}_\mu^{-1}$ with standard BPX: \begin{align} \mathbf{\tilde{D}}_\mu^{-1} = \sum_{k=0}^L \mathbf{R}_k \left(h_k^{-1} \mathbf{M}_k + \mu \mathbf{M}_k \right)^{-1} \mathbf{R}_k^T \end{align}

• $L$ nested meshes with mesh width $h_k$.
• $\mathbf{M}_k$ is the standard mass matrix on level $k$.
• $\mathbf{R}_k$ is the embedding matrix from level $k$ onto the finest level $L$.
• $\mathbf{\tilde{D}}_\mu^{-1}$ is spectrally equivalent, so the symmetric scaling proof stays valid.

\begin{align} u(t) = \sin(\pi t) \in C^\infty(\overline{I}), \quad \mu=10 \end{align}

• Convergence order $p_t + 0.5$ in the energy norm and $p_t + 1$ in the $L^2$ norm.

\begin{align} u(t) = \sqrt{t} \in H^{1-\varepsilon}_{0,}(I), \quad \mu=10 \end{align}

• Convergence order $0.5$ in the energy and $L^2$ norms.

For $u \in \mathcal{Q}_s = H^{1/2+\varepsilon}_{0,}(I; L^2(\Omega)) \cap L^2(I; H^1_0(\Omega))$, we consider the variational problem

\begin{align} b_h(u,v) = \langle \partial_t u, (I + \mathcal{H}_T) v \rangle_Q + \langle \nabla_x u, (I + \mathcal{H}_T) \nabla_x v \rangle_Q = \langle f, (I + \mathcal{H}_T) v \rangle_Q \quad \forall v \in \mathcal{Q}. \end{align}

• $\mathcal{Q} = H^{1/2}_{0,}(I; L^2(\Omega)) \cap L^2(I; H^1_0(\Omega))$ is the Lions-Magenes space.
• Ellipticity with respect to $\mathcal{Q}$ and boundedness with respect to $\mathcal{Q}_s$ and $\mathcal{Q}$.
• Quasi-best approximation can be used for error estimates.
• Theory works for simplicial and tensor product meshes.
• On tensor product meshes ($\mathcal{Q}_h = S_{h_t; 0,}^{p_t} \otimes X_h$) $\qquad \mathbf{B} \mathbf{u} = \mathbf{f}, \qquad \mathbf{B} = \mathbf{A}_t \otimes \mathbf{M}_x + \mathbf{M}_t \otimes \mathbf{K}_x$
  FFT techniques as well as BPX preconditioning remain applicable.

$I=(0,1)$, $\Omega=(0,1)^2$, uniform meshes with $h_t=h_x$.

• For further levels, uniform refinement is used.

\begin{align} u(t,x,y) = \sin(\pi t) \sin(\pi x) \sin(\pi y) \in C^\infty(\overline{Q}). \end{align} Order $p$ in the energy norm and $p+1$ in the $L^2$ norm.

\begin{align} u(t,x,y) = \sqrt{t} \sin(\pi x) \sin(\pi y) \in H^{1-\varepsilon}(I, C^\infty(\overline{\Omega})). \end{align} Only $p=1$ is meaningful here, consistent with the ODE case.
Reduced convergence can be cured using parabolic scaling.