Continuous time integration for changing type systems

TL;DR

This paper introduces a continuous time integration method for systems with evolving types, proving optimal convergence and demonstrating its effectiveness through numerical examples.

Contribution

It develops a variational time integration approach using continuous Galerkin Petrov methods for changing type systems, with proven optimal error convergence.

Findings

Proven optimal-order convergence of the method.

Numerical examples demonstrating the method's effectiveness.

Theoretical analysis of error bounds.

Abstract

We consider variational time integration using continuous Galerkin Petrov methods applied to evolutionary systems of changing type. We prove optimal-order convergence of the error in a cGP-like norm and conclude the paper with some numerical examples and conclusions.

Tables

Table 1. Table 1: Errors and rates for example ( 5.1 )

	cGP-method				dG-method
$M=2N$	${\Vert |U_{ref}-U_h^{\tau }|\Vert }_{\rho }$		${\Vert U_{ref}-U_h^{\tau }\Vert }_{\rho ,[0,T]}$		${\Vert U_{ref}-U_h^{\tau }\Vert }_{\rho ,[0,T]}$
$k=2$, $r=1$
256	2.120e-02		8.890e-04		1.808e-04
512	5.746e-03	1.88	3.136e-04	1.50	7.751e-05	1.22
1024	1.787e-03	1.68	1.380e-04	1.18	3.496e-05	1.15
2048	7.036e-04	1.35	6.739e-05	1.03	1.580e-05	1.15
$k=3$, $r=2$
256	8.806e-04	1.05	1.187e-04		6.058e-05
512	4.163e-04	1.08	5.489e-05	1.11	2.642e-05	1.20
1024	1.906e-04	1.13	2.492e-05	1.14	1.137e-05	1.22
2048	8.581e-05	1.15	1.114e-05	1.16	4.669e-06	1.28

Table 2. Table 2: Errors ‖ U r ​ e ​ f − U h τ ‖ ρ , [ 0 , T ] subscript norm subscript 𝑈 𝑟 𝑒 𝑓 subscript superscript 𝑈 𝜏 ℎ 𝜌 0 𝑇 \|{U_{ref}-U^{\tau}_{h}}\|_{\rho,[0,T]} and rates for example ( 5.2 )

	cGP-method				dG-method
$M=2N$	${\Vert |U_{ref}-U_h^{\tau }|\Vert }_{\rho }$		${\Vert U_{ref}-U_h^{\tau }\Vert }_{\rho ,[0,T]}$		${\Vert U_{ref}-U_h^{\tau }\Vert }_{\rho ,[0,T]}$
$k=2$, $r=1$
16	3.989e-02		1.961e-02		7.821e-03
32	1.972e-02	1.02	9.199e-03	1.09	3.018e-03	1.37
64	9.435e-03	1.06	3.751e-03	1.29	8.813e-04	1.78
96	5.603e-03	1.29	1.324e-03	1.50	2.920e-04	1.59
$k=3$, $r=2$
16	1.041e-02		5.499e-03		2.790e-03
32	3.689e-03	1.50	1.435e-03	1.94	6.385e-04	2.13
64	1.248e-03	1.56	4.430e-04	1.70	2.248e-04	1.51

Equations

$${\partial }_t{U}_1+\mathrm{div}(U_2)$$

$$U_1+\mathrm{div}(U_2)$$

$$M_0=\left(\begin{array}{cc}[c]{\chi }_{{\Omega }_{\mathrm{hyp}}\cup {\Omega }_{\mathrm{par}}}& 0\\ 0& {\chi }_{{\Omega }_{\mathrm{hyp}}}\end{array}\right),M_1=\left(\begin{array}{cc}[c]{\chi }_{{\Omega }_{\mathrm{ell}}}& 0\\ 0& {\chi }_{{\Omega }_{\mathrm{par}}\cup {\Omega }_{\mathrm{ell}}}\end{array}\right)\text{and}A=\left(\begin{array}{cc}[c]0& \mathrm{div}\\ {\mathrm{grad}}^{\circ }& 0\end{array}\right),$$

$$({\partial }_t{M}_0+M_1+A)U=F.$$

$$M_{0}U(0^{+})=M_0{U}_0.$$

$$H_{\rho }(\mathbb{R};\mathbf{H}):=\left\{f:\mathbb{R}\to \mathbf{H}:f\text{mbox}meas.,{\int }_{\mathbb{R}}\Vert f(t){\Vert }_{\mathbf{H}}^2\exp (-2\rho t)\mathrm{d}t<\infty \right\}.$$

$$⟨f,g{⟩}_{\rho }:={\int }_{\mathbb{R}}⟨f(t),g(t)⟩\exp (-2\rho t)\mathrm{d}t$$

$$⟨(\rho M_0+M_1)x,x⟩\ge \gamma ⟨x,x⟩=\gamma \Vert x{\Vert }_{\mathbf{H}}^2,$$

$$F{|}_{{\mathbb{R}}_{\ge 0}}\text{ is continuous and }F(t)=0,t<0,$$

$$U_0\in \text{dom}(A)$$

$$(M_1+A)U_0=F(0^{+})$$

$$U(0^{+})=U_0.$$

$$({\partial }_t{M}_0+M_1+A)U=F+{\delta }_0{M}_0{U}_0$$

$$({\partial }_t{M}_0+M_1+A)(U-H_0{U}_0)$$

$$\mathbf{H}=L^2(\Omega )\otimes (L^2(\Omega ){)}^n\text{and}D(A)=H_0^1(\Omega )\otimes H_{\mathrm{div}}(\Omega ).$$

$${\mathcal{U}}_h^{\tau }$$

$${\mathcal{V}}_h^{\tau }$$

$$V_1$$

$$V_2$$

$$V_1\otimes V_2\subset D(A)\subset \mathbf{H}.$$

$$R{T}_{k-1}(\sigma )=({\mathcal{P}}_{k-1}(\sigma ){)}^n+\mathbf{x}{\mathcal{P}}_{k-1}(\sigma )\subset {\mathcal{P}}_k(\sigma {)}^n.$$

$$⟨f,g{⟩}_{\rho ,m}:={\int }_{I_m}⟨f(t),g(t)⟩\exp (-2\rho t)\mathrm{d}t$$

$$B_m(U_h^{\tau },V_h^{\tau })$$

$$U_h^{\tau }(0)$$

$$⟨U-{\Pi }_h^{\tau }U,W_h^{\tau }{⟩}_{\rho }=0,\text{for all }{W}_h^{\tau }\in {\mathcal{V}}_h^{\tau },$$

$${\left|\left|\left|U_h^{\tau }\right|\right|\right|}_{\rho }^2:=\frac{1}{2}\Vert M_0^{1/2}{U}_h^{\tau }(T){\Vert }_{\mathbf{H}}^2{\mathrm{e}}^{-2\rho T}+\Vert N{U}_h^{\tau }(0){\Vert }_{\mathbf{H}}^2+\gamma \Vert {\Pi }_h^{\tau }{U}_h^{\tau }{\Vert }_{\rho }^2.$$

$$N{U}_h^{\tau }(0)=0$$

$$R{U}_h^{\tau }(T)=0$$

$$U_h^{\tau }=R{U}_h^{\tau }+N{U}_h^{\tau }=0$$

$$\frac{1}{2}\Vert M_0^{1/2}{U}_h^{\tau }(T){\Vert }_{\mathbf{H}}^2{\mathrm{e}}^{-2\rho T}+\gamma \Vert {\Pi }_h^{\tau }{U}_h^{\tau }{\Vert }_{\rho }^2\le \sum _{m=1}^M{B}_m(U_h^{\tau },{\Pi }_h^{\tau }{U}_h^{\tau })+\frac{1}{2}\Vert M_0^{1/2}\mathcal{I}{U}_0{\Vert }_{\mathbf{H}}^2.$$

$$B_m(U_h^{\tau },{\Pi }_h^{\tau }{U}_h^{\tau })=⟨{\partial }_t{M}_0{U}_h^{\tau },U_h^{\tau }{⟩}_{\rho ,m}+⟨M_1{\Pi }_h^{\tau }{U}_h^{\tau },{\Pi }_h^{\tau }{U}_h^{\tau }{⟩}_{\rho ,m},$$

$$\langle\partial_{t}M_{0}U^{\tau}_{h},U^{\tau}_{h}\rangle_{\rho,m}=\rho\langle M_{0}U^{\tau}_{h},U^{\tau}_{h}\rangle_{\rho,m}+\frac{1}{2}\|{M_{0}^{1/2}U^{\tau}_{h}(t)}\|_{\mathbf{H}}^{2}\mathrm{e}^{-2\rho t}\big{|}_{t_{m-1}}^{t_{m}}.$$

$$⟨M_0{U}_h^{\tau },U_h^{\tau }{⟩}_{\rho ,m}$$

Full text

AMS subject classification (2010): 65J08, 65J10, 65M12, 65M60

Key words: evolutionary equations, changing type system, continuous Galerkin Petrov, space-time approach

1 Introduction

Let us start with an example, where the type of the problem changes over the spacial domain and has homogeneous Dirichlet boundary conditions. For this purpose let $n\in\{1,2,3\}$ be the spatial dimension and $\Omega\subset\mathbb{R}^{n}$ be bounded and partitioned into measurable, disjoint sets $\Omega_{\mathrm{ell}},\Omega_{\mathrm{par}}$ and $\Omega_{\mathrm{hyp}}$. In $\Omega_{\mathrm{hyp}}$ a hyperbolic wave equation is given for $U=(U_{1},U_{2})$

[TABLE]

Each of above equations is also known in their derived second order formulation for $U_{1}$, namely $(\partial_{t}^{2}-\operatorname{\Delta})U_{1}=\partial_{t}F_{1}-\operatorname{div}F_{2}$ for the wave equation, $(\partial_{t}-\operatorname{\Delta})U_{1}=F_{1}-\operatorname{div}F_{2}$ for the heat equation and $(1-\operatorname{\Delta})U_{1}=F_{1}-\operatorname{div}F_{2}$ for the reaction-diffusion equation.

Denoting by $\chi_{D}$ the characteristic function of a domain $D\subset\Omega$ and defining the linear operators

[TABLE]

where ∘ denotes the homogeneous Dirichlet boundary conditions w.r.t. $\Omega$, we can write above equations in a condensed way

[TABLE]

Now we are left with the question, under which conditions above problem (1.1) has a unique solution.

In the following we assume $U_{0}$ in $D(A)$. Besides that condition we can draw a condition on the operators from a much more general theory. 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 [7] that this form is an evolutionary problem, given by (1.1), where $\partial_{t}$ stands for the derivative with respect to time, $M_{0}:\mathbf{H}\to\mathbf{H}$ and $M_{1}:\mathbf{H}\to\mathbf{H}$ are bounded linear selfadjoint operators on some Hilbert space $\mathbf{H}$, $A:D(A)\subset\mathbf{H}\to\mathbf{H}$ is an unbounded skew-selfadjoint operator on $\mathbf{H}$ and $F$ is a given source term.

We are interested in a unique solution $U$ of above equation. For this purpose let $\rho>0$ and define the weighted $L^{2}$-function space

[TABLE]

The space $H_{\rho}(\mathbb{R};\mathbf{H})$ is a Hilbert space endowed with the natural inner product given by

[TABLE]

for all $f,g\in H_{\rho}(\mathbb{R};\mathbf{H})$, where $\left\langle f(t),g(t)\right\rangle$ is the inner product of $\mathbf{H}$ and $\|{\cdot}\|_{\mathbf{H}}$ its associated norm. We obtain a norm by setting $\|{f}\|_{\rho}^{2}:=\langle f,f\rangle_{\rho}$. The associated weighted $H^{k}$-function spaces are denoted by $H_{\rho}^{k}(\mathbb{R};\mathbf{H})$ for $k\in\mathbb{N}$. Now from [7, Thm. (solution theory)] it follows: If there exists a $\rho_{0}>0$ and a $\gamma>0$ such that for all $\rho\geq\rho_{0}$ and $x\in\mathbf{H}$

[TABLE]

then for all right hand sides $F\in H_{\rho}(\mathbb{R},\mathbf{H})$ exists a unique solution $U\in H_{\rho}(\mathbb{R},\mathbf{H})$. Furthermore, by above condition $\left\langle M_{0}x,x\right\rangle\geq 0$ follows and there exists a root $M_{0}^{1/2}$ of $M_{0}$. Note that the theory presented in [7] deals with vanishing initial conditions at $t\to-\infty$.

Corollary 1.1**.**

Under the conditions (1.2) and

[TABLE]

problem (1.1) has a unique solution $U$ with

[TABLE]

Proof.

Problem (1.1) given as a problem on $\mathbb{R}$ reads

[TABLE]

where the initial condition $M_{0}U(0^{+})=M_{0}U_{0}$ is included via the delta distribution $\delta_{0}$ at $t=0$ on the right-hand side.

Let $H_{0}$ denote the Heaviside function with the jump at $t=0$. We obtain for $U-H_{0}U_{0}$ the evolutionary problem

[TABLE]

By (1.3) we have $\tilde{F}(t)=0,\,t<0$, $\tilde{F}(0)=0$ and $\tilde{F}$ is continuous. Now [10] yields that the problem for $U-H_{0}U_{0}$ has a unique solution in $H_{\rho}^{1}(\mathbb{R},\mathbf{H})$. Thus $U$ is a unique solution of (1.1) and $U(0^{+})=U_{0}$. ∎

In the following we assume conditions (1.2) and (1.3) to be fulfilled. Then $U_{0}$ is an initial data on the whole $\Omega$, explicitly also in the elliptic and parabolic regime. But due to the compatibility condition (1.3) it cannot be chosen independently of $F$.

In [6] the class of changing type problems 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 [11, 4, 8, 1, 2, 3].

Note that, like in [6], we deal in this paper with problems that have a changing type over the given domain and could be rewritten into second order form as shown above. But then transmission conditions would need to be stated that are embedded automatically into the first order formulation. This is a very useful feature of the general approach and it allows to combine models from different parts of physics into one well-posed problem. We want to emphasise that the time discretisation presented and analysed in this paper holds for all problems of above general class of first order problems, only the spatial discretisation has to be adapted to the operator $A$.

For our problem and operator $A$ the Hilbert space $\mathbf{H}$ and $D(A)$ can now be specified to

[TABLE]

Remark 1.2**.**

*The solution theory demands $A$ to be skew-selfadjoint which in turn restricts the choice of boundary data. Some simple choices are homogeneous Dirichlet boundary conditions on the first component, encoded by $\operatorname{grad}^{\circ}$ in above operator $A$, homogeneous Neumann boundary conditions on the second component, encoded by $\operatorname{div}^{\circ}$ or periodic boundary conditions on both components, encoded by $\operatorname{grad}^{\#}$ and $\operatorname{div}^{\#}$.

Inhomogeneous conditions can always be transformed into homogeneous ones by a substitution changing the right hand side of the problem.*

The paper is organised as follows. The precise formulation of the method considered is stated in Section 2 while Section 3 deals with the existence of discrete solutions. In Section 4 we present error estimates and finally Section 5 gives some numerical examples and conclusions.

2 Numerical method

The discrete variational form of (1.1) 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\}$. Furthermore 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}}$, i.e. each of these subdomains is a union of mesh cells, and let $h$ be the maximal diameter of the cells of $\Omega_{h}$. Furthermore, let $r,\,k\geq 1$ denote polynomial degrees.

Then the piecewise polynomial function spaces for the trial and test functions resp. are given by

[TABLE]

where the spatial spaces are

[TABLE]

and therefore

[TABLE]

Here $\mathcal{P}_{k}(\sigma)$ is the space of polynomials of degree up to $k$ on the cell $\sigma$ of $\Omega_{h}$ and $RT_{k-1}(\sigma)$ is the Raviart-Thomas-space, defined by

[TABLE]

Note that we retain the regularity in space of the trial functions also for the test functions in order to define a Galerkin method in space. Furthermore, if the mesh consists of quadrilateral or hexahedral cells, in above definitions and statements the polynomial space $\mathcal{P}_{k}(\sigma)$ can be replaced by a mapped $\mathcal{Q}_{k}$-space, including all polynomials of total degree $k$ over a reference element mapped onto $\sigma$. If the mesh is a combination of both types of cells, a combination of spaces also works with a suitable mapping ensuring the continuities.

Let us localise in addition the scalar product in $H_{\rho}(\mathbb{R},\mathbf{H})$ to the time intervals $I_{m}$ by

[TABLE]

and the norm $\|{f}\|_{\rho,m}^{2}:=\langle f,f\rangle_{\rho,m}$. Then the variational formulation using the continuous Galerkin-Petrov method reads:

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\}$

[TABLE]

is the initial value. Here $\mathcal{I}=(\mathcal{I}_{1},\mathcal{I}_{2})$ denotes the spatial interpolation operator, where $\mathcal{I}_{1}:H_{\rho}([0,T],H^{1}(\Omega))\to H_{\rho}([0,T],V_{1})$ is locally the Scott–Zhang interpolant on each cell $\sigma$, see [9] for a precise definition, and $\mathcal{I}_{2}:H_{\rho}((0,t),H(\operatorname{div},\Omega)\cap(L^{s}(\Omega))^{n})\to H_{\rho}([0,T],V_{2})$ with $s>2$ is the standard interpolator defined via moments, see [5]. Note that it is appropriate to include the full initial conditions into the discrete problem, see Corollary 1.1.

3 Existence of discrete solution

Let us start by defining $\Pi_{h}^{\tau}$ as the orthogonal $L^{2}$-projection w.r.t. $\langle\cdot,\cdot\rangle_{\rho}$ into the test space $\mathcal{V}^{\tau}_{h}$, i.e.

[TABLE]

$R$ and $N$ as the projectors onto the range and nullspace of $M_{0}$, resp, and

[TABLE]

Lemma 3.1**.**

The seminorm $\left|\!\!\;\left|\!\!\;\left|{U^{\tau}_{h}}\right|\!\!\;\right|\!\!\;\right|_{\rho}$ is a norm on $\mathcal{U}^{\tau}_{h}$.

Proof.

With $\mathcal{U}^{\tau}_{h}$ being finite we only have to show that $\left|\!\!\;\left|\!\!\;\left|{U^{\tau}_{h}}\right|\!\!\;\right|\!\!\;\right|_{\rho}=0$ implies $U^{\tau}_{h}=0$. Thus, let us assume $\left|\!\!\;\left|\!\!\;\left|{U^{\tau}_{h}}\right|\!\!\;\right|\!\!\;\right|_{\rho}=0$. Then it follows immediately $\Pi_{h}^{\tau}U^{\tau}_{h}=0$ and due to continuity, one degree of freedom is left for $U^{\tau}_{h}$. On each time interval $U^{\tau}_{h}$ is a multiple of a weighted Legendre polynomial that is orthogonal to $V^{\tau}_{h}$ w.r.t. $\langle\cdot,\cdot\rangle_{\rho}$. From $\|{NU^{\tau}_{h}(0)}\|_{\mathbf{H}}=0$ we conclude

[TABLE]

and therefore $NU^{\tau}_{h}=0$, because the Legendre polynomial is not zero at the left boundary. From $\|{M_{0}^{1/2}U^{\tau}_{h}(T)}\|_{\mathbf{H}}=0$ we have similarly

[TABLE]

and therefore $RU^{\tau}_{h}=0$, because the Legendre polynomial is not zero at the right boundary. With

[TABLE]

we have the assertion. ∎

Lemma 3.2**.**

It holds

[TABLE]

Proof.

Let us consider any interval $I_{m}$. Then it holds

[TABLE]

where the skew-symmetry of $A$ and the definition of $\Pi_{h}^{\tau}$ was used. For the first term we apply integration by parts and obtain due to the exponential weight

[TABLE]

By the $L^{2}$-orthogonality (3.1) it follows

[TABLE]

and therefore

[TABLE]

With the general existence assumption $\rho M_{0}+M_{1}\geq\gamma$ and $M_{0}\geq 0$ we obtain

[TABLE]

Summing over the intervals the statement follows. ∎

It follows

[TABLE]

and therefore

[TABLE]

This shows unique existence and continuous dependence on $f$ and $U_{0}$ of the discrete solution $U^{\tau}_{h}$.

4 Error-estimation

Let us start by stating interpolation error estimates.

Interpolation in time

Let $P_{r}:H_{\rho}^{1}([0,T],\mathbf{H})\to H^{1}_{\rho}([0,T],\mathbf{H})$, where P_{r}u\big{|}_{I_{m}}\in\mathcal{P}_{r}(I_{m},\mathbf{H}) for all $m\in\{1,\dots,M\}$, be the interpolation operator fulfilling locally for all $m$ and $v\in H_{\rho}^{1}([0,T],\mathbf{H})$

[TABLE]

Although we have weighted norms and scalar products the standard interpolation error estimates

[TABLE]

holds for $v\in H^{r+1}_{\rho}([0,T],\mathbf{H})$, where here and further on $C>0$ denotes a generic constant and $\tau:=\max\{\tau_{m}\}$.

Interpolation in space

As previously stated we use $\mathcal{I}=(\mathcal{I}_{1},\mathcal{I}_{2})$ as spatial interpolation operator, where the first component $\mathcal{I}_{1}:H_{\rho}([0,T],H^{1}(\Omega))\to H_{\rho}([0,T],V_{1})$ is the Scott–Zhang interpolant, and the second component $\mathcal{I}_{2}:H_{\rho}((0,t),H(\operatorname{div},\Omega)\cap(L^{\sigma}(\Omega))^{n})\to H_{\rho}([0,T],V_{2})$ with $\sigma>2$ is the standard Raviart-Thomas interpolator. Here it holds for all $v\in H_{0}^{1}(\Omega)\cap H^{s}(\Omega)$, see [9],

[TABLE]

where $1\leq s\leq k+1$, $\|{v}\|_{s}$ denotes the $H^{s}(\Omega)$-norm, and for all $q\in H^{s}(\Omega)$ such that $\operatorname{div}q\in H^{s}(\Omega)$, see [5]

[TABLE]

where $1\leq s\leq k$.

Error analysis

Note that we do have for all $V^{\tau}_{h}\in\mathcal{V}^{\tau}_{h}$ the Galerkin orthogonality

[TABLE]

for the solution $U\in H_{\rho}^{1}([0,T],\mathbf{H})$ of (1.1) and $U^{\tau}_{h}\in\mathcal{U}^{\tau}_{h}$ of (2.1). We now want to estimate the error $U-U^{\tau}_{h}$ and decompose it into $U-U^{\tau}_{h}=\eta+\xi$, where

[TABLE]

Note that with (2.1b) it follows

[TABLE]

Lemma 4.1**.**

It holds for any $m\in\{1,\dots,M\}$ and $V^{\tau}_{h}\in\mathcal{V}^{\tau}_{h}$

[TABLE]

Proof.

Using the Galerkin orthogonality (4.3) we obtain the error equality

[TABLE]

Using integration by parts and the properties of $P_{r}$ we obtain for all $w\in\mathcal{V}^{\tau}_{h}$ and $v\in H_{\rho}^{1}([0,T],\mathbf{H})$

[TABLE]

Thus we get the error equation

[TABLE]

from which (4.4) follows by a Cauchy-Schwarz inequality. ∎

From the error equation (4.5) and the stability estimate (3.3) we obtain

[TABLE]

by substituting $U^{\tau}_{h}:=\xi$ and $f:=-(2\rho M_{0}+M_{1}+A)\eta$, and noting $\xi(0)=0$.

In order to simplify the representation of the main result, let us abbreviate

[TABLE]

and

[TABLE]

Theorem 4.2**.**

We assume for the solution $U$ of (1.1) the regularity

[TABLE]

as well as

[TABLE]

Then we have for the error of the numerical solution $U^{\tau}_{h}$ of (2.1)

[TABLE]

Proof.

By the decomposition of the norm and the error we have to estimate

[TABLE]

Using above interpolation error estimates we obtain

[TABLE]

For the remaining two terms we apply (4.6) and obtain

[TABLE]

and similarly for $\|{M_{0}^{1/2}\xi(T)}\|_{H}$. Combining these results proves the error estimate. ∎

Remark 4.3**.**

We assumed in Theorem 4.2 slightly higher regularity assumptions on $U$ than actually needed. Instead of assuming $U\in H^{1}_{\rho}([0,T],\mathbf{H}^{k})$ for the point evaluation at $t=T$ the weaker assumption $U\in W^{0,\infty}_{\rho}([0,T],\mathbf{H}^{k})$ suffices. But in order to prove that claim from conditions on the right-hand side the easiest way is by proving above regularity and using the Sobolev-embedding.

Remark 4.4**.**

*In this section we presented an error analysis for the fully discrete problem of the changing type system. At the same time it is true for all operators $M_{0}$ and $M_{1}$ fulfilling Assumption (1.2). The analysis can also easily be adapted to general evolutionary problems having a different spatial operator $A$ by defining suitable discrete spatial function spaces and corresponding interpolation operators, and providing sufficient interpolation error estimates. *

Theorem 4.5**.**

In the case of $M_{0}>0$, e.g. a purely hyperbolic problem, we can also give a convergence result in the weighted $L^{2}$-type $\|{\cdot}\|_{\rho}$-norm. Under the same conditions as in Theorem 4.2 we have

[TABLE]

Proof.

For this result we need

• •

a local norm equivalence for all $W_{h}^{\tau}\in\mathcal{U}_{h}^{\tau}$

[TABLE]

with a constant $C_{1}$ independent of $\tau_{m}$ and $W^{\tau}_{h}$, that holds true because $\Pi_{h}^{\tau}W^{\tau}_{h}-W^{\tau}_{h}$ is a multiple of a weighted Legendre polynomial of degree $r$, $t_{m}$ is not a zero of it and the scaling w.r.t. $\tau_{m}$ of the two terms is the same,

• •

a local estimation of the discrete error $\xi$ with a localisation of the norms to the interval $[0,t_{m}]$ instead of $[0,T]$

[TABLE]

that follows by the same lines as in Thm. 4.2,

• •

a Sobolev embedding for $t_{n}<t_{m}$ and $U\in H^{1}_{\rho}([t_{n},t_{m}],\mathbf{H})$

[TABLE]

with a constant $C_{inv}$ independent of $U$, $t_{n}$ and $t_{m}$.

Then it follows

[TABLE]

where the Sobolev embedding for $\|{U(T)}\|_{\mathbf{H}^{k}}\mathrm{e}^{-\rho T}$ uses the whole interval $[0,T]$ and only $[t_{m-1},t_{m}]$ for $\|{U(t_{m})}\|_{\mathbf{H}^{k}}\mathrm{e}^{-\rho t_{m}}$, as well as $\tau_{m}\leq 1$. Together with the interpolation error bound

[TABLE]

the claim follows. ∎

5 Numerical examples

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. That both examples have initial values zero is not a restriction. We look into the convergence behaviour also w.r.t. the weighted $L^{2}$-norm $\|{\cdot}\|_{\rho}$ in addition to the $\left|\!\!\;\left|\!\!\;\left|{\cdot}\right|\!\!\;\right|\!\!\;\right|_{\rho}$-norm in order to compare the results with those of the discontinuous Galerkin method from [6]. In the finite discrete setting both norms are equivalent.

All computations were done in the finite-element framework $\mathbb{SOFE}$222github.com/SOFE-Developers/SOFE.

5.1 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_{\mathrm{hyp}}=\left[-\pi,0\right]$, and $\Omega_{\mathrm{ell}}=\left[0,\pi\right]$. As final time we set $T=4\pi$. The problem is now given by

[TABLE]

with homogeneous Dirichlet-conditions for the first component of $U:\mathbb{R}\times\mathbb{R}\to\mathbb{R}\times\mathbb{R}$, the initial condition $U_{0}=0$ and a right-hand side $F(t,x)=(f(t,x),g(t,x))\cdot\chi_{\geq 0}(t)$, where $\chi_{\geq 0}(t)$ is the characteristic function of the non-negative time line and

[TABLE]

Thus, $F$ is continuous on $\mathbb{R}$ and it holds $F(t)=0$ for $t\leq 0$. Therefore, the solution theory of [7] gives the existence of a unique solution $U$ that is continuous in time. Figure 1

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\times 2048$ mesh with polynomial degrees $k=4$ and $r=3$. Table 1

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 only a convergence rate between 1 and 2 in both norms. Increasing the polynomial degree reduces the error, but does not improve the rate much. A reason for this behaviour could be that $U$ is not smooth enough for the error estimates to hold due having jumping coefficients in space and a non-differentiable right hand side. Unfortunately the exact solution to this problem and thus its precise regularity is unknown.

For comparison we also computed approximations with the discontinuous Galerkin method from [6] that uses globally discontinuous piecewise polynomials of degree $r$ in time and the same approximation in space as the method described in this paper. 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.

5.2 1+2d example

As second example we consider the last example of [6]. 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

[TABLE]

where

[TABLE]

Figure 2

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 2

shows the results. Similarly to the previous example we do not achieve the optimal convergence order for both methods. Here the data and the right-hand side have jumps along interior lines which reduces the maximum regularity of the solution. Again the discontinuous Galerkin method has smaller errors.

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 in the $\left|\!\!\;\left|\!\!\;\left|{\cdot}\right|\!\!\;\right|\!\!\;\right|_{\rho}$-norm. 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 enough, 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.

Acknowledgement

The author would like to thank R. Picard for many helpful discussions on the existence of solution for evolutionary problems and the anonymous reviewer for the helpful comments. Also thanks goes to L. Ludwig for providing $\mathbb{SOFE}$.