\section{Introduction}
  Most of the classical linear partial differential equations arising in mathematical physics
  can be written in a common operator form. It has been shown in \cite{Picard} that this form
  is an evolutionary problem, given by
  \begin{equation}\label{eq:evo}
    (\partial_tM_0+M_1+A)U=F,
  \end{equation}
  where $\partial_t$ stands for the derivative with respect to time, 
  $M_0$ and $M_1$ are bounded linear operators on some Hilbert space $H$,
  $A$ is an unbounded skew-selfadjoint operator on $H$
  and $F$ is a given source term. 
  
  We are interested in a unique solution $U$ of above equation. The theory presented in \cite{Picard}
  deals vanishing initial conditions at $t\to-\infty$, but any non-zero initial condition
  can be incorporated as data into $F$.
  For this purpose let $\rho> 0$ and define the weighted $L^2$-function space 
  \[
    H_\rho(\mathbb{R};H):=
      \left\{ 
          f:\mathbb{R}\to H\,:\, f \mbox{ meas.}, 
          \int_\mathbb{R} \norm{f(t)}{H}^2 \exp(-2\rho t) \,\mathrm{d}t<\infty
      \right\}.
  \]
  The space $H_\rho(\mathbb{R};H)$ is a Hilbert space endowed with the natural inner 
  product given by 
  \[
    \scpr{f,g}\,:=
    \int_{\mathbb{R}} \scp{f(t),g(t)} \exp(-2\rho t) \,\mathrm{d}t
  \]
  for all $f,g\in H_\rho(\mathbb{R};H)$, where $\scp{f(t),g(t)}$ is the inner product of $H$. 
  The associated weighted $H^k$-function spaces are denoted by $H_\rho^k(\mathbb{R};H)$ for $k\in \mathbb{N}$.
  Now the solution theory of \cite{Picard} states: If there exists a $\rho_0>0$ and a $\gamma>0$
  such that
  \[
    \rho M_0+M_1\geq \gamma\quad\text{for all }\rho\geq\rho_0,
  \]
  then for all right hand sides $F\in H_{\rho}(\mathbb{R},H)$ exists a unique solution $U\in H_{\rho}(\mathbb{R},H)$.
  
  In \cite{FrTW16} this class of problems was investigated numerically using
  a discontinuous Galerkin approach for the discretisation in time. Here we want to apply 
  a continuous approach, namely the continuous Galerkin-Petrov method \cite{Winther81,AM89,Schieweck10,AM15}.
  Like in \cite{FrTW16} we will consider a special subclass of the general class of evolutionary problems,
  although the approach of the time-discretisation is suitable in the general case.
  
  For the precise class of problems let $n\in\{1,2,3\}$ be the spatial dimension and 
  $\Omega\subset \mathbb{R}^n$ be partitioned into measurable, disjoint sets 
  $\Omega_{\mathrm{ell}},\Omega_{\mathrm{par}}$ and $\Omega_{\mathrm{hyp}}$. Denoting by $\id{D}$
  the characteristic function of a domain $D\subset\Omega$ we define the linear operators 
  \[
    M_0=\begin{pmatrix}[c]
            \id{\Omega_{\mathrm{hyp}}\cup\Omega_{\mathrm{par}}} & 0\\
            0 &\id{\Omega_{\mathrm{hyp}}}
        \end{pmatrix}\quad\text{and}\quad
    M_1=\begin{pmatrix}[c]
          \id{\Omega_{\mathrm{ell}}} & 0 \\
          0 & \id{\Omega_{\mathrm{par}}\cup\Omega_{\mathrm{ell}}}
        \end{pmatrix}.     
  \]
  Additionally let the unbounded skew-selfadjoint operator A be given by  
  \[
    A=\begin{pmatrix}[c]
          0 & \operatorname{div}\\
          \operatorname{grad}^\circ &0
       \end{pmatrix}.    
  \]
  Here $^\circ$ denotes the homogeneous Dirichlet boundary conditions. The Hilbert space $H$ can now be 
  specified to 
  \[
    H(\Omega)=H_0^1(\Omega)\otimes H_{\operatorname{div}}(\Omega).
  \]
        
  The resulting evolutionary problem is now of mixed type. More precisely, 
  on $\Omega_{\mathrm{ell}}$ we get an equation of elliptic type, 
  on $\Omega_{\mathrm{par}}$ the equations becomes parabolic while 
  on $\Omega_{\mathrm{hyp}}$ the problem is hyperbolic.

  The solution theory of \cite{Picard} does not only give for a final time $T>0$ and $F\in H_\rho([0,T],H)$ 
  the existence of a unique solution $U\in H_\rho([0,T],H)$.
  It holds in addition $M_0U\in H_\rho^1([0,T],H)$ and therefore, by the Sobolev embedding theorem, $M_0U$ is continuous in time.
  If the right-hand side is smoother, say $F\in H_\rho^1([0,T],H)$, and $F$ 
  and $(M_1+A)U_0$ are in the range of $M_0$, then the continuity in time of the solution 
  is everywhere and we have, see \cite{TW19}
  \[
    U\in H_\rho^1([0,T],H)\subset C([0,T],H).
  \]
  In the following we want to utilise this additional regularity.
  Section~\ref{sec:method} gives the precise formulation of the method considered.
  Section~\ref{sec:exist} then deals with the existence of discrete solutions, 
  Section~\ref{sec:error} presents error estimates and finally
  Section~\ref{sec:numerics} gives some numerical examples and conclusions.
  
\section{Numerical method}\label{sec:method}
  Let us start wit a semi-discretisation in time.
  The semi-discrete variational form of \eqref{eq:evo} uses a decomposition of $[0,T]$ into $M$
  disjoint intervals $I_m=(t_{m-1},t_m]$ of length $\tau_m=t_m-t_{m-1}$ for $m\in\{1,\dots,M\}$
  and the piecewise polynomial-in-time spaces
  \begin{align*}
    \mathcal{U}^\tau&:=\{u\in H^1_\rho(\mathbb{R},H):u\big|_{I_m}\in\mathcal{P}_{r}(I_m,H),m\in\{1,\dots,M\}\},\\
    \mathcal{V}^\tau&:=\{v\in H_\rho(\mathbb{R},H):v\big|_{I_m}\in\mathcal{P}_{r-1}(I_m,H),m\in\{1,\dots,M\}\}
  \end{align*}
  for the trial and the test functions, where $\mathcal{P}_r(I_m,H)$ is the space of polynomials of degree up to $r$
  on $I_m$ with values in $H$. Let us localise the scalar product in $H_\rho(\mathbb{R},H)$
  to the time intervals $I_m$ by
  \[
    \scprm{f,g}\,:=
    \int_{I_m} \scp{f(t),g(t)} \exp(-2\rho t) \,\mathrm{d}t.
  \]  
  Then the variational formulation using the continuous Galerkin-Petrov method reads:\\
  Find $U^\tau\in\mathcal{U}^\tau$ such that for all $V^\tau\in\mathcal{V}^\tau$ and $m\in\{1,\dots,M\}$
  \begin{align}
    \scprm{(\partial_t M_0+M_1+A)U^\tau,V^\tau}&=\scprm{f,V^\tau}=\scprm{\Pi_{r-1}f,V^\tau},\label{eq:evosemidisc}\\
    U^\tau(t_{m-1}^+)&=U^\tau(t_{m-1}^-),\notag
  \end{align}
  where $\Pi_{r-1}$ denotes the weighted $L^2$-projection into $\mathcal{V}^\tau$ and $U^\tau(0^-)=U_0$ is 
  the initial value. 
  
  In order to compute an approximation to the solution of problem \eqref{eq:evo}, we need 
  also to discretise the spatial domain and the function space $H$. Let $\Omega$ be discretised 
  into $\Omega_h$ by a regular simplicial mesh that resolves the sets $\Omega_{\mathrm{ell}},\,
  \Omega_{\mathrm{par}}$ and $\Omega_{\mathrm{hyp}}$, and let $h$ be the maximal diameter 
  of the cells of $\Omega_h$. Furthermore, let $k\geq 0$ denote a polynomial degree.
  Then the discrete spaces are given by
  \begin{align*}
    \mathcal{U}^\tau_h 
    \!&:=\! \left\{
          u_h\in\mathcal{U}^\tau\!:
          u_h|_{I_m}\!\in\mathcal{P}_r(I_m,V_1\otimes V_2),
          m\!\in\!\{1,\dots,M\}
       \right\},\\
    \mathcal{V}^\tau_h 
    \!&:=\! \left\{
          v_h\in\mathcal{V}^\tau\!:
          v_h|_{I_m}\!\in\mathcal{P}_{r-1}(I_m,V_1\otimes V_2),
          m\!\in\!\{1,\dots,M\}
       \right\},
  \end{align*}
  where the spatial spaces are
  \begin{align*}
    V_1&:=\left\{v\in H_0^1(\Omega):\,v|_\sigma\in\mathcal{P}_k(\sigma)\,\forall \sigma\in\Omega_h\right\},\\
    V_2&:=\left\{w\in H_{\operatorname{div}}(\Omega):\,w|_\sigma\in RT_{k-1}(\sigma)\,\forall \sigma\in\Omega_h\right\}.
  \end{align*}
  Here $\mathcal{P}_k(\sigma)$ is the space of polynomials of degree up to $k$ on 
  the cell $\sigma$ and $RT_{k-1}(\sigma)$ is the Raviart-Thomas-space, defined by
  \[
    RT_{k-1}(\sigma)=(\mathcal{P}_{k-1}(\sigma))^n+\mathbf{x}\mathcal{P}_{k-1}(\sigma)\subset\mathcal{P}_k(\sigma)^n.
  \]
  The fully discrete method then reads similarly to \eqref{eq:evosemidisc}:\\
  Find $U^\tau_h\in\mathcal{U}^\tau_h$ such that for all $V^\tau_h\in\mathcal{V}^\tau_h$ and $m\in\{1,\dots,M\}$
  \begin{align}
    B_m(U^\tau_h,V^\tau_h):=
    \scprm{(\partial_t M_0+M_1+A)U^\tau_h,V^\tau_h}&=\scprm{f,V^\tau_h},\label{eq:evodisc}\\
    U^\tau_h(t_{m-1}^+)&=U^\tau_h(t_{m-1}^-),\notag
  \end{align}
  where and $U^\tau_h(0^-)=U_0$ is the initial value.
  
\section{Existence of discrete solution}\label{sec:exist}
  Let us start by choosing $V^\tau_h=\Pi_{r-1}U^\tau_h$, the $L^2$-projection of $U^\tau_h$ 
  into the test space $\mathcal{V}^\tau_h$ and denote by 
  \[
    \norm{f}{\rho,[0,t_\ell]}^2:=\sum_{m=1}^\ell\scprm{f,f}
  \]
  a localised version of $\norm{f}{\rho}^2=\scpr{f,f}$.

  \begin{lemma}\label{lem:exist}
    It holds for all $\ell\geq 1$
    \[
      \scp{M_0 U^\tau_h(t_\ell),U^\tau_h(t_\ell)}_H\mathrm{e}^{-2\rho t_\ell}+
      \gamma\norm{\Pi_{r-1}U^\tau_h}{\rho,[0,t_\ell]}^2
      \leq\sum_{k=1}^\ell B_k(U^\tau_h,\Pi_{r-1}U^\tau_h)+\scp{M_0 U_0,U_0}_H.
    \]
  \end{lemma}
  \begin{proof}
    Let us consider any interval $I_m$. Then it holds
    \[
      B_m(U^\tau_h,\Pi_{r-1}U^\tau_h)
        = \scprm{\partial_t M_0 U^\tau_h,U^\tau_h}+\scprm{M_1 \Pi_{r-1}U^\tau_h,\Pi_{r-1}U^\tau_h},
    \]
    where the skew-symmetry of $A$ and the definition of $\Pi_{r-1}$ was used.
    
    For the first term we apply integration by parts and obtain due to the exponential weight
    \[
      \scprm{\partial_t M_0 U^\tau_h,U^\tau_h}
        =\rho\scprm{M_0 U^\tau_h,U^\tau_h}+\scp{M_0 U^\tau_h(t),U^\tau_h(t)}_H\mathrm{e}^{-2\rho t}\big|_{t_{m-1}}^{t_m}.
    \]
    By the $L^2$-orthogonality it follows
    \begin{align*}
      \scprm{M_0 U^\tau_h,U^\tau_h}
        &= \scprm{M_0 (U^\tau_h-\Pi_{r-1}U^\tau_h),U^\tau_h-\Pi_{r-1}U^\tau_h}+\scprm{M_0 \Pi_{r-1}U^\tau_h,\Pi_{r-1}U^\tau_h}\\
        &\geq \scprm{M_0 \Pi_{r-1}U^\tau_h,\Pi_{r-1}U^\tau_h}
    \end{align*}
    and therefore
    \[
      \scprm{\partial_t M_0 U^\tau_h,U^\tau_h}
        \geq \scprm{\rho M_0 \Pi_{r-1}U^\tau_h,\Pi_{r-1}U^\tau_h}+\scp{M_0 U^\tau_h(t),U^\tau_h(t)}_H\mathrm{e}^{-2\rho t}\big|_{t_{m-1}}^{t_m}.
    \]
    With the general existence assumption $\rho M_0+M_1\geq \gamma$ 
    and $M_0\geq 0$ we obtain
    \begin{gather}
      B_m(U^\tau_h,\Pi_{r-1}U^\tau_h)
        \geq \gamma\norm{\Pi_{r-1}U^\tau_h}{\rho,m}^2+
              \scp{M_0 U^\tau_h(t),U^\tau_h(t)}_H\mathrm{e}^{-2\rho t}\big|_{t_{m-1}}^{t_{m}}.\label{eq:lowerbound}
    \end{gather}
    Summing over the intervals the statement follows.
  \end{proof}
  With Lemma~\ref{lem:exist} we have control over the $L^2$-projection of the solution $U^\tau_h$ and of the final value of $\norm{M_0^{1/2}U^\tau_h(t)}{H}$. 
  In order to obtain control of the full solution $U$ we consider two approaches.

  \textbf{1) $\mathbf{M_0}$ has a trivial null space}\\
  In this case we employ the norm equivalence in discrete spaces
  \[
    \norm{U^\tau_h}{\rho,m}^2\sim\norm{\Pi_{r-1}U^\tau_h}{\rho,m}^2+\tau_m\scp{M_0 U^\tau_h(t_m),U^\tau_h(t_m)}_H\mathrm{e}^{-2\rho t_m}
  \]
  with constants independent of $\tau_m$ and $U^\tau_h$ and we obtain 
  \begin{align*}
    \norm{U^\tau_h}{\rho,[0,T]}^2
      &\leq C\left[\norm{\Pi_{r-1}U^\tau_h}{\rho,[0,T]}^2
        +\sum_{m=1}^M\tau_m\scp{M_0U^\tau_h(t_m),U^\tau_h(t_m)}_H\mathrm{e}^{-2\rho t_m}\right]\\
      &\leq C\left(\frac{1}{\gamma}+T\right)\left[\sum_{m=1}^MB_m(U^\tau_h,\Pi_{r-1}U^\tau_h)+\scp{M_0 U_0,U_0}_H\right]\\
      &\leq C\left(\frac{1}{\gamma}+T\right)\left[\scp{\Pi_{r-1}f,U^\tau_h}_{\rho,[0,T]}+\scp{M_0 U_0,U_0}_H\right].
  \end{align*}
  Applying a Young inequality we have
  \begin{align*}
    \norm{U^\tau_h}{\rho,[0,T]}^2
      &\leq C^2\left(\frac{1}{\gamma}+T\right)^2\norm{f}{\rho,[0,T]}^2+
                2C\left(\frac{1}{\gamma}+T\right)\scp{M_0 U_0,U_0}_H
  \end{align*}
  Therefore, in the case of $M_0$ with only trivial null-space the linear system is solvable 
  and we also have continuous dependence of the solution $U^\tau_h$ on the data. 

  \textbf{2) $\mathbf{M_0}$ may have a non-trivial null space}\\
  Here we can apply another norm-equivalence to show the solvability of the discrete systems. 
  Using $\rho M_0+M_1\geq\gamma>0$ we have the norm equivalence 
  \[
    \norm{U^\tau_h}{\rho,m}^2\sim\norm{\Pi_{r-1}U^\tau_h}{\rho,m}^2+\tau_m\scp{(\rho M_0+M_1) U^\tau_h(t_{m-1}),U^\tau_h(t_{m-1})}_H\mathrm{e}^{-2\rho t_{m-1}}
  \]
  with constants independent of $\tau_m$ and $U^\tau_h$.
  Using a similar approach as before we obtain
  \begin{align*}
    \norm{U^\tau_h}{\rho,[0,T]}^2
      &\leq C\Bigg[\frac{1}{\gamma}\left(\sum_{m=1}^MB_m(U^\tau_h,\Pi_{r-1}U^\tau_h)+\scp{M_0 U_0,U_0}_H\right)\\
      &\qquad      +\sum_{m=1}^{M-1}\rho\tau_{m+1}\left(\sum_{n=1}^m B_n(U^\tau_h,\Pi_{r-1}U^\tau_h)+\scp{M_0 U_0,U_0}_H\right)+\rho\tau_1\scp{M_0U_0,U_0}_H\\
      &\qquad      +\sum_{m=1}^M\tau_m\scp{M_1 U^\tau_h(t_{m-1}),U^\tau_h(t_{m-1})}_H\mathrm{e}^{-2\rho t_{m-1}}\Bigg]\\
      &\leq C\Bigg[\frac{1}{\gamma}\left(\scp{\Pi_{r-1}f,U^\tau_h}_{\rho,[0,T]}+\scp{M_0 U_0,U_0}_H\right)\\
      &\qquad      +\sum_{m=1}^{M-1}\rho\tau_{m+1}\left(\scp{\Pi_{r-1}f,U^\tau_h}_{\rho,[0,t_m]}+\scp{M_0 U_0,U_0}_H\right)+\rho\tau_1\scp{M_0U_0,U_0}_H\\
      &\qquad      +\sum_{m=1}^M\tau_m\scp{M_1 U^\tau_h(t_{m-1}),U^\tau_h(t_{m-1})}_H\mathrm{e}^{-2\rho t_{m-1}}\Bigg]\\
      &\leq C\Bigg[\left(\frac{1}{\gamma}+\rho T\right)\norm{\Pi_{r-1}f}{\rho,[0,T]}\norm{U^\tau_h}{\rho,[0,T]}+\left(\frac{1}{\gamma}+\rho (T+\tau_1)\right)\scp{M_0 U_0,U_0}_H\\
      &\qquad      +T\max_{m\in\{1,\dots,M\}}\scp{M_1 U^\tau_h(t_{m-1}),U^\tau_h(t_{m-1})}_H\mathrm{e}^{-2\rho t_{m-1}}\Bigg].
  \end{align*}
  With a Young inequality follows
  \begin{multline*}
    \norm{U^\tau_h}{\rho,[0,T]}^2
      \leq C^2\left(\frac{1}{\gamma}+\rho T\right)^2\norm{f}{\rho,[0,T]}^2+2C\left(\frac{1}{\gamma}+\rho (T+\tau_1)\right)\scp{M_0 U_0,U_0}_H\\
      \qquad      +2CT\max_{m\in\{1,\dots,M\}}\scp{M_1 U(t_{m-1}),U(t_{m-1})}_H\mathrm{e}^{-2\rho t_{m-1}}.
  \end{multline*}
  Therefore, also in the general case the linear system is solvable 
  and the solution on each time interval depends continuously on the data and the initial value of the time interval. 
  
\section{Error-estimation}\label{sec:error}
  \subsection{Semi discretisation in time}
    Note that we do have the Galerkin orthogonality
    \begin{gather}\label{eq:Galorth}
      B_m(U-U^\tau,V)=0
    \end{gather}
    for the solution $U\in H_\rho^1(\mathbb{R},H)$ of \eqref{eq:evo} for $F\in H_\rho^1(\mathbb{R},H)$ 
    and $U^\tau\in \mathcal{U}^\tau$ of \eqref{eq:evosemidisc}. We now want to estimate the error $U-U^\tau$.    
    Let $P_r:H_\rho^1([0,T],H)\to\mathcal{U}^\tau$ be the interpolation operator fulfilling locally 
    for all $m$ and $v\in H_\rho^1([0,T],H)$
    \begin{align*}
      (P_rv-v)(t_{m-1})&=0,\\
      (P_rv-v)(t_{m})&=0,\\
      \scprm{P_ru-u,w}&=0\quad\forall w\in\mathcal{P}_{r-2}(I_m,H).
    \end{align*}
    We decompose the error into $u-U=\eta+\xi$, where
    \[
      \eta=u-P_ru,\quad
      \xi=P_ru-U\in\mathcal{U}^\tau.
    \]
    \begin{lemma}
      It holds for any $m\in\{1,\dots,M\}$
      \begin{gather}\label{eq:errineq}
        \scprm{(\partial_t M_0+M_1+A)\xi,V}
          \leq\left(\norm{(2\rho M_0+M_1)\eta}{\rho,m}+
                    \norm{A\eta}{\rho,m}\right)\norm{V}{\rho,m}.
      \end{gather}      
    \end{lemma}
    \begin{proof}
      Using the Galerkin orthogonality \eqref{eq:Galorth} it follows
      the error equality
      \[
        \scprm{(\partial_t M_0+M_1+A)\xi,V}
          =-\scprm{(\partial_t M_0+M_1+A)\eta,V}.
      \]
      Then we have by integration by parts and the properties of $P_r$ for all $w\in\mathcal{V}^\tau$ and $v\in H_\rho^1([0,T],H)$
      \begin{align*}
        \scprm{\partial_tM_0(v-P_r v),w}
          &= 2\rho\scprm{M_0(v-P_r v),w}-\scprm{v-P_r v,\partial_t M_0w}+\scp{v-P_r v,w}_H\mathrm{e}^{-2\rho t}\big|_{t_{m-1}}^{t_m}\\
          &= 2\rho\scprm{M_0(v-P_r v),w}.
      \end{align*}
      Thus we get the error equation
      \[
        \scprm{(\partial_t M_0+M_1+A)\xi,V}
          =-\scprm{(2\rho M_0+M_1+A)\eta,V}
      \]
      from which \eqref{eq:errineq} follows by a Cauchy-Schwarz inequality.
    \end{proof}
    In the following estimates we will bound the terms of interest by a norm of the interpolation error given by
    \[
      \norm{\eta}{*,[0,t_m]}:=\norm{(2\rho M_0+M_1)\eta}{\rho,[0,t_m]}+
                              \norm{A\eta}{\rho,[0,t_m]}.
    \]    
    \begin{lemma}
      It holds for any $m\in\{1,\dots,M\}$
      \begin{align*}
        \gamma\norm{\Pi_{r-1}\xi}{\rho,[0,t_m]}^2+
        \scp{M_0\xi(t_m),\xi(t_m)}_H\mathrm{e}^{-2\rho t_m}
        &\leq\frac{2}{\gamma}\norm{\eta}{*,[0,t_m]}^2      
      \end{align*}
    \end{lemma}
    \begin{proof}
      Let $V^\tau=\Pi_{r-1}\xi$ be the test function in \eqref{eq:errineq}. 
      From \eqref{eq:lowerbound} and 
      \eqref{eq:errineq}
      we have
      \begin{align*}
        \gamma\norm{\Pi_{r-1}\xi}{\rho,m}^2+\frac{1}{2}\scp{M_0\xi(t),\xi(t)}_H\mathrm{e}^{-2\rho t}\big|_{t_{m-1}}^{t_m}
         &\leq
        \scprm{(\partial_t M_0+M_1+A)\xi,\Pi_{r-1}\xi}\\
        &\leq\left(\norm{(2\rho M_0+M_1)\eta}{\rho,m}+
                  \norm{A\eta}{\rho,m}\right)\norm{\Pi_{r-1}\xi}{\rho,m}.        
      \end{align*}
      Now Young's-inequality ($ab\leq \frac{1}{2\gamma}a^2+\frac{\gamma}{2}b^2$) leads to
      \[
        \gamma\norm{\Pi_{r-1}\xi}{\rho,m}^2+\scp{M_0\xi(t),\xi(t)}_H\mathrm{e}^{-2\rho t}\big|_{t_{m-1}}^{t_m}
        \leq\frac{1}{\gamma}\left(\norm{(2\rho M_0+M_1)\eta}{\rho,m}+
                                  \norm{A\eta}{\rho,m}\right)^2
      \]
      which gives upon summation the statement because $M_0\xi(0)=0$.       
    \end{proof}
    Thus we have control over the weighted $L^2$-projection of the error 
    and its final value.
    \begin{corollary}
      In the case of $M_0$ having only a trivial null-space we have immediate control of the full norm
      due to
      \[
        \norm{\xi}{\rho,m}^2\leq C\left(\norm{\Pi_{r-1}\xi}{\rho,m}^2+\tau_m\scp{M_0\xi(t_m),\xi(t_m)}_H\mathrm{e}^{-2\rho t_m}\right)
      \]
      for some $C>0$ implying
      \[
        \norm{\xi}{\rho,[0,T]}^2
          \leq C \left(\frac{2}{\gamma^2}+\frac{2T}{\gamma}\right)\norm{\eta}{*,[0,T]}^2.
      \]
      In the case of $M_0$ having a non-trivial null-space the norm equivalence considered in 
      the existence proof of Section~\ref{sec:exist} do only provide a bound exponentially increasing 
      in $1/\min\tau_m$.
    \end{corollary}
    
  \subsection{Full discretisation}
    Let $\tau:=\max\limits_{m\in\{1,\dots,M\}}\tau_m$ be the maximal time-step length.
    Similarly to \cite[Section 4]{FrTW16} we can now proceed and estimate the error of
    the full discretisation. Here the only difference is the time-discretisation method.
    \begin{theorem}\label{thm:cts}
      We assume for the solution $U$ of \eqref{eq:evo} the 
      regularity     
      \[
        U\in H_\rho^{1}(\mathbb{R};H^k(\Omega)\times H^k(\Omega))\cap 
        H_\rho^{r+3}(\mathbb{R};L^2(\Omega)\times L^2(\Omega)) 
      \]
      as well as 
      \[
        AU\in H_\rho(\mathbb{R}; H^k(\Omega)\times H^k(\Omega)).
      \]
      Then we have for the error of the numerical solution $U^\tau_h$ by 
      \eqref{eq:evodisc}
      \[
        \norm{U-\Pi_{r-1}U_h^\tau}{\rho,[0,T]}^2+
        \mathrm{e}^{-2\rho T}\scp{M_0(U-U_h^\tau)(T),(U-U_h^\tau)(T)}
        \leq C (\tau^{2(r+1)} + T h^{2k})
      \]
      and if $M_0$ has only a trivial null space
      \[
        \norm{U-U_h^\tau}{\rho,[0,T]}^2
        \leq C (\tau^{2(r+1)} + T h^{2k}).
      \]
    \end{theorem}
    \begin{proof}
      This is a direct consequence of the results in \cite[Section 4]{FrTW16} 
      combined with the estimates of the previous section.
    \end{proof}

\section{Numerical examples}\label{sec:numerics}
  We consider two examples with unknown solutions. Simulations with known smooth solutions 
  were also made and the theoretical orders were observed. The two following examples 
  show a more realistic behaviour in the case of changing type systems.
  All computations were done in the finite-element framework $\mathbb{SOFE}$, see 
  \texttt{github.com/SOFE-Developers/SOFE}.
  \subsection{1+1d example}
    Let us consider as first example one spatial dimension and combine a hyperbolic and an elliptic region.
    To be more precise, let $\Omega=\left[-\pi,\pi\right]$, $\Omega_{hyp}= \left[-\pi,0\right]$, 
    and $\Omega_{ell}=\left[0,\pi\right]$. As final time we set $T=4\pi$. The problem is now given
    by
    \begin{gather}\label{ex:1d}
      \left[
        \partial_t
        \begin{pmatrix}[c]
        \id{\Omega_{hyp}} & 0\\
        0 & \id{\Omega_{hyp}}
        \end{pmatrix}
        +\begin{pmatrix}[c]
          \id{\Omega_{ell}} & 0\\
          0 & \id{\Omega_{ell}}
        \end{pmatrix}
        +\begin{pmatrix}[c]
          0 & \partial_x\\
          \mathring{\partial_x} & 0
        \end{pmatrix}
      \right]U
        =F
    \end{gather}
    with homogeneous Dirichlet-conditions for the first component of $U:\mathbb{R}\times\mathbb{R}\to\mathbb{R}^2$,
    the initial condition $U_0=0$
    and a right-hand side $F(t,x)=(f(t,x),g(t,x))\cdot\id{\geq 0}(t)$, where $\id{\geq 0}(t)$ 
    is the characteristic function of the non-negative time line and
    \begin{align*}
      f(t,x) &= \frac{1}{5}\sin(3t)+\min\{t,\pi\}\cos(3x),\\
      g(t,x) &= \sin(t)\left(1-\frac{x^2}{\pi^2}\right).
    \end{align*}
    Thus, $F$ is continuous and it holds $F(t)=0$ for $t\leq 0$. 
    Therefore, the solution theory of \cite{Picard} gives the existence of a 
    unique solution $U$ that is continuous in time. Figure~\ref{fig:sol1d}
    \begin{figure}
     \begin{center}
        \includegraphics[width=0.35\textwidth]{ex1da}
        \includegraphics[width=0.35\textwidth]{ex1db}
      
      
     \end{center}
      \caption{Solution of problem \eqref{ex:1d}, first component (left) and second component (right)\label{fig:sol1d}}
    \end{figure}
    shows plots of the components of the solution in the domain. Note that 
    the first component has a kink along $x=0$ -- it is continuous but not 
    differentiable in $x$. As mesh we use an equidistant mesh of $N$ cells in $\Omega$
    and $M$ cells in $[0,T]$. In order to calculate the errors we use a reference 
    solution $U_{ref}$ instead of the unknown solution $U$. The reference solution is computed 
    on an $4096\times2048$ mesh with polynomial degrees $k=4$ and $r=3$.
    Table~\ref{tab:ex1d}
    \begin{table}
      \caption{Errors $\norm{U_{ref}-U^\tau_h}{\rho,[0,T]}$ and rates for example \eqref{ex:1d}\label{tab:ex1d}}
      \begin{center}
        \begin{tabular}{rllll}
          \toprule
          $M=2N$
          &\multicolumn{2}{c}{cGP-method}
          &\multicolumn{2}{c}{dG-method}\\
          \midrule
          \multicolumn{5}{c}{$k=2$, $r=1$}\\
          \midrule
            256 & 8.890e-04 &      & 1.808e-04 &     \\
            512 & 3.136e-04 & 1.50 & 7.751e-05 & 1.22\\
           1024 & 1.380e-04 & 1.18 & 3.496e-05 & 1.15\\
           2048 & 6.739e-05 & 1.03 & 1.580e-05 & 1.15\\
          \midrule
          \multicolumn{5}{c}{$k=3$, $r=2$}\\
          \midrule
            256 & 1.187e-04 &      & 6.058e-05 &     \\
            512 & 5.489e-05 & 1.11 & 2.642e-05 & 1.20\\
           1024 & 2.492e-05 & 1.14 & 1.137e-05 & 1.22\\
           2048 & 1.114e-05 & 1.16 & 4.669e-06 & 1.28\\
           \bottomrule
        \end{tabular}        
      \end{center}
    \end{table}
    shows the results for different values of $M$ and $N$ and polynomial degrees $k$ and $r$.
    We coupled $k=r+1$ as the theory gives for smooth $U$ the convergence order $\min\{k,r+1\}$
    if $N$ and $M$ are proportional. We observe for the continuous Galerkin-Petrov method 
    in the first column only a convergence rate between 1 and 2. Increasing the polynomial
    degree reduces the error, but does not influence the rate much. A reason for this behaviour 
    could be that $U$ is not smooth enough for the error estimates to hold.\\
    For comparison we also computed approximations with the discontinuous Galerkin method 
    from \cite{FrTW16}. The errors given in the remaining columns show a similar behaviour 
    with convergence rates between 1 and 2. Nevertheless, the errors are smaller 
    for the discontinuous approach.
    
  \subsection{1+2d example}
    As second example we consider the last example of \cite{FrTW16}.
    Let $T=5.2$, $\Omega=(0,1)^2\subset\mathbb{R}^2$, 
    $\Omega_{\mathrm{hyp}}=\left( \frac{1}{4},\frac{3}{4}\right)^2$ and 
    $\Omega_{\mathrm{ell}}=\Omega\setminus\bar\Omega_{\mathrm{hyp}}$
    The problem is given by
    \begin{gather}\label{ex:2d}
      \left[
      \partial_t
      \begin{pmatrix}
        \id{\Omega_{\mathrm{hyp}}} & 0\\
        0 & \id{\Omega_{\mathrm{hyp}}}
      \end{pmatrix}
      +\begin{pmatrix}
        \id{\Omega_{\mathrm{ell}}} & 0\\
        0 & \id{\Omega_{\mathrm{ell}}}
      \end{pmatrix}
      +\begin{pmatrix}
        0 & \operatorname{div} \\
        \operatorname{grad}^\circ & 0
      \end{pmatrix}
      \right]U
      =\begin{pmatrix}
        f\\
        0
      \end{pmatrix},
    \end{gather}
    where
    \[
      f(t,\mathbf{x})=2\sin(\pi t)\cdot\id{\mathbb{R}_{<1/2}\times\mathbb{R}}(\mathbf{x}).
    \]
    Figure~\ref{fig:sol2d} 
 
    \begin{figure}[tb]
      \begin{center}
        \includegraphics[width=0.25\textwidth]{ex2d1.pdf}\quad
        \includegraphics[width=0.25\textwidth]{ex2d2.pdf}\quad
        \includegraphics[width=0.25\textwidth]{ex2d3.pdf}\\
        \includegraphics[width=0.25\textwidth]{ex2d4.pdf}\quad
        \includegraphics[width=0.25\textwidth]{ex2d5.pdf}\quad
        \includegraphics[width=0.25\textwidth]{ex2d6.pdf}
      \end{center}
      \caption{First component of $U$ at times $t=5k/16$ for $k\in\{1,\dots, 6\}$ (top left 
               to bottom right) of problem \eqref{ex:2d}}\label{fig:sol2d}
    \end{figure}
   
    shows some snapshots of the first component of the solution $U:\mathbb{R}\times\mathbb{R}^2\to\mathbb{R}\times\mathbb{R}^2$,
    approximated by a numerical simulation. Again we use equidistant meshes
    with $N$ cells in each dimension of space and $M$ cells in $[0,T]$.
    As reference solution $U_{ref}$ replacing the unknown exact solution we use
    an approximation calculated with $M=192,\,N=96,\,k=3,\,r=2$ and $M=128,\,N=64,\,k=4,\,r=3$, resp.
    Table~\ref{tab:ex2d}
    \begin{table}
      \caption{Errors $\norm{U_{ref}-U^\tau_h}{\rho,[0,T]}$ and rates for example \eqref{ex:2d}\label{tab:ex2d}}
      \begin{center}   
        \begin{tabular}{rllll}
          \toprule
          $M=2N$
          &\multicolumn{2}{c}{cGP-method}
          &\multicolumn{2}{c}{dG-method}\\
          \midrule
          \multicolumn{5}{c}{$k=2$, $r=1$}\\
          \midrule
          16 & 1.961e-02 &      & 7.821e-03 &     \\
          32 & 9.199e-03 & 1.09 & 3.018e-03 & 1.37\\
          64 & 3.751e-03 & 1.29 & 8.813e-04 & 1.78\\
          96 & 1.324e-03 & 1.50 & 2.920e-04 & 1.59\\
          \midrule
          \multicolumn{5}{c}{$k=3$, $r=2$}\\
          \midrule
          16 & 5.499e-03 &      & 2.790e-03 &     \\
          32 & 1.435e-03 & 1.94 & 6.385e-04 & 2.13\\
          64 & 4.430e-04 & 1.70 & 2.248e-04 & 1.51\\
          \bottomrule
        \end{tabular}
      \end{center}    
    \end{table}
    shows the results. Similarly to the previous example we do not achieve the optimal 
    convergence order for both methods. And again the discontinuous Galerkin method has smaller errors.
    
\section*{Conclusions}
  The continuous solution of an evolutionary system with continuous right hand side
  can be approximated by several methods. Here we investigated the continuous Galerkin-Petrov
  method, that has optimal convergence order for smooth solutions. We have proved this
  for the $L^2$-projection of the error into the test space. For operators $M_0$ having 
  only a trivial null space it also follows for the full error. The benefit of the continuous 
  method compared to the discontinuous Galerkin method is the continuity that implies a 
  non-dissipative behaviour. In our examples with unknown solutions, that are probably not smooth,
  the discontinuous Galerkin method is slightly better. Furthermore, these examples show that 
  an increase of the polynomial degree in space over 2 and in time over 1 gives no huge benefit. 
  This is different for smooth solutions -- here both methods achieve the theoretical high convergence orders.
  

\bibliographystyle{plain}
