Numerical study of Galerkin-collocation approximation in time for the wave equation

TL;DR

This paper introduces a higher order Galerkin-collocation method for time discretization of wave equations, combining accuracy with computational efficiency, and validates its performance through numerical experiments.

Contribution

It establishes a novel connection between Galerkin and collocation methods for wave equations, enabling high-order accuracy with reduced computational costs.

Findings

Achieves high-order temporal accuracy in wave simulations.

Reduces computational complexity compared to traditional methods.

Demonstrates improved efficiency and accuracy through numerical tests.

Abstract

The elucidation of many physical problems in science and engineering is subject to the accurate numerical modelling of complex wave propagation phenomena. Over the last decades, high-order numerical approximation for partial differential equations has become a well-established tool. Here we propose and study numerically the implicit approximation in time of wave equations by a Galerkin--collocation approach that relies on a higher order space-time finite element approach. The conceptual basis is the establishment of a direct connection between the Galerkin method for the time discretization and the classical collocation methods, with the perspective of achieving the accuracy of the former with reduced computational costs provided by the latter in terms of less complex linear algebraic systems. For the fully discrete solution, higher order regularity in time is further ensured which can be advantageous in the discretization of multi-physics systems. The accuracy and efficiency of the variational collocation approach is carefully studied by numerical experiments.

Figures

Tables

Table 1. Table 1: Calculated errors for GCC ( 3 ) 1 {}^{1}(3) with solution ( 39 ).

$\tau$	$h$	${\Vert e^u\Vert }_{L^{\infty }\left(L^2\right)}$	EOC	${\Vert e^v\Vert }_{L^{\infty }\left(L^2\right)}$	EOC	${\Vert \left|E\right|\Vert }_{L^{\infty }}$	EOC
${\tau }_0/2^0$	$h_0$	2.318e-04	–	1.543e-03	–	1.574e-03	–
${\tau }_0/2^1$	$h_0$	1.541e-05	3.91	9.694e-05	3.99	1.004e-04	3.97
${\tau }_0/2^2$	$h_0$	9.825e-07	3.97	6.260e-06	3.95	6.478e-06	3.95
${\tau }_0/2^3$	$h_0$	6.185e-08	3.99	3.946e-07	3.99	4.082e-07	3.99
${\tau }_0/2^4$	$h_0$	3.873e-09	4.00	2.472e-08	4.00	2.557e-08	4.00
${\tau }_0/2^5$	$h_0$	2.422e-10	4.00	1.548e-09	4.00	1.609e-09	3.99
$\tau$	$h$	${\Vert e^u\Vert }_{L^2\left(L^2\right)}$	EOC	${\Vert e^v\Vert }_{L^2\left(L^2\right)}$	EOC	${\Vert \left|E\right|\Vert }_{L^2}$	EOC
${\tau }_0/2^0$	$h_0$	1.634e-04	–	1.232e-03	–	1.441e-03	–
${\tau }_0/2^1$	$h_0$	1.070e-05	3.93	7.864e-05	3.97	9.269e-05	3.96
${\tau }_0/2^2$	$h_0$	6.765e-07	3.98	4.943e-06	3.99	5.836e-06	3.99
${\tau }_0/2^3$	$h_0$	4.240e-08	4.00	3.094e-07	4.00	3.654e-07	4.00
${\tau }_0/2^4$	$h_0$	2.652e-09	4.00	1.934e-08	4.00	2.285e-08	4.00
${\tau }_0/2^5$	$h_0$	1.659e-10	4.00	1.212e-09	4.00	1.433e-09	3.99

Table 2. Table 2: Calculated errors for GCC ( 3 ) 1 {}^{1}(3) with solution ( 40 ).

$\tau$	$h$	${\Vert e^u\Vert }_{L^{\infty }\left(L^2\right)}$	EOC	${\Vert e^v\Vert }_{L^{\infty }\left(L^2\right)}$	EOC	${\Vert \left|E\right|\Vert }_{L^{\infty }}$	EOC
${\tau }_0/2^0$	$h_0/2^0$	3.486e-03	–	3.602e-02	–	5.013e-02	–
${\tau }_0/2^1$	$h_0/2^1$	2.329e-04	3.90	2.392e-03	3.90	3.338e-03	3.92
${\tau }_0/2^2$	$h_0/2^2$	1.483e-05	3.97	1.527e-04	3.97	2.128e-04	3.98
${\tau }_0/2^3$	$h_0/2^3$	9.320e-07	3.99	9.609e-06	3.99	1.338e-05	3.99
${\tau }_0/2^4$	$h_0/2^4$	5.837e-08	4.00	6.022e-07	4.00	8.383e-08	4.00
${\tau }_0/2^5$	$h_0/2^5$	3.649e-09	4.00	3.767e-08	4.00	5.243e-08	4.00
$\tau$	$h$	${\Vert e^u\Vert }_{L^2\left(L^2\right)}$	EOC	${\Vert e^v\Vert }_{L^2\left(L^2\right)}$	EOC	${\Vert \left|E\right|\Vert }_{L^2}$	EOC
${\tau }_0/2^0$	$h_0/2^0$	2.700e-03	–	2.568e-02	–	3.458e-02	–
${\tau }_0/2^1$	$h_0/2^1$	1.771e-04	3.93	1.689e-03	3.93	2.278e-03	3.92
${\tau }_0/2^2$	$h_0/2^2$	1.120e-05	3.98	1.070e-04	3.98	1.444e-04	3.98
${\tau }_0/2^3$	$h_0/2^3$	7.020e-07	4.00	6.713e-06	3.99	9.061e-06	3.99
${\tau }_0/2^4$	$h_0/2^4$	4.391e-08	4.00	4.199e-07	4.00	5.669e-07	4.00
${\tau }_0/2^5$	$h_0/2^5$	2.744e-09	4.00	2.624e-08	4.00	3.543e-08	4.00

Table 3. Table 3: Runtime (wall clock time) for GCC ( 3 ) 1 {}^{1}(3) (method C 1 superscript 𝐶 1 C^{1} ) and cGP(2) (method C 0 superscript 𝐶 0 C^{0} ) for different time step sizes and solvers of section   4.2.1 ( r 1 subscript 𝑟 1 r_{1} ) and 4.2.2 ( r 2 subscript 𝑟 2 r_{2} ) .

DoF (space)	cores	method	${\tau }_n$	$r_1$[h]	$r_2$[h]
$3.2\text{×}{10}^6$	224	$C^0$	$0.25\times {\tau }_0$	$219.3$	-
		$C^0$	${\tau }_0$	$40.0$	-
		$C^1$	${\tau }_0$	$46.6$	$25.3$
		$C^1$	$2\times {\tau }_0$	$33.1$	$19.4$
		$C^1$	$25\times {\tau }_0$	$4.5$	$2.3$
		$C^1$	$35\times {\tau }_0$	$3.7$	$2.2$
		$C^1$	$50\times {\tau }_0$	$2.9$	$1.6$
		$C^1$	$100\times {\tau }_0$	$1.7$	$0.9$
		$C^1$	$200\times {\tau }_0$	$1.1$	$0.7$
$4.2\text{×}{10}^6$	336	$C^1$	$50\times {\tau }_0$	$3.3$	$1.7$

Table 4. Table 4: Calculated errors for GCC ( 5 ) 2 {}^{2}(5) with solution ( 39 ).

$\tau$	$h$	${\Vert e^u\Vert }_{L^{\infty }\left(L^2\right)}$	EOC	${\Vert e^v\Vert }_{L^{\infty }\left(L^2\right)}$	EOC	${\Vert \left|E\right|\Vert }_{L^{\infty }}$	EOC
${\tau }_0/2^0$	$h_0$	8.748e-06	–	4.355e-05	–	4.985e-05	–
${\tau }_0/2^1$	$h_0$	1.370e-07	6.00	7.404e-07	5.88	8.043e-07	5.95
${\tau }_0/2^2$	$h_0$	2.165e-09	5.98	1.202e-08	5.95	1.266e-08	5.99
${\tau }_0/2^3$	$h_0$	3.388e-11	6.00	1.883e-10	6.00	1.980e-10	6.00
${\tau }_0/2^4$	$h_0$	5.301e-13	6.00	2.940e-12	6.00	3.093e-12	6.00
$\tau$	$h$	${\Vert e^u\Vert }_{L^2\left(L^2\right)}$	EOC	${\Vert e^v\Vert }_{L^2\left(L^2\right)}$	EOC	${\Vert \left|E\right|\Vert }_{L^2}$	EOC
${\tau }_0/2^0$	$h_0$	4.022e-06	–	2.996e-05	–	3.502e-05	–
${\tau }_0/2^1$	$h_0$	6.353e-08	5.98	4.808e-07	5.96	5.599e-07	5.97
${\tau }_0/2^2$	$h_0$	9.957e-10	6.00	7.565e-09	5.99	8.800e-09	5.99
${\tau }_0/2^3$	$h_0$	1.557e-11	6.00	1.184e-10	6.00	1.377e-10	6.00
${\tau }_0/2^4$	$h_0$	2.431e-13	6.00	1.849e-12	6.00	2.151e-12	6.00

Equations

$$\begin{array}[]{r@{\;}c@{\;}ll}\partial_{t}^{2}u-c^{2}\Delta{u}&=f\,,&\in\Omega\times I\,,\\[4.30554pt] u(0)=u_{0}\,,\quad\partial_{t}u(0)&=v_{0}\,,&\text{in}\;\Omega\,,\\[4.30554pt] u=g^{u}\,,\;\;\text{on}\;\partial\Omega_{D}\times I\,,\quad\partial_{n}u&=0\,,&\text{on}\;\partial\Omega_{N}\times I\,.\\ \end{array}$$

$${\partial }_{t}u-v=0,{\partial }_{t}v-c^2\Delta u=f.$$

$$u=u^0+u^D\text{and}v=v^0+v^D.$$

$$u^0(0)=u_0-u^D(0),v^0(0)=v_0-v^D(0)$$

$$\displaystyle\int_{I}\bigl{<}\partial_{t}u^{0},\phi\bigr{>}_{\Omega}-\bigl{<}v^{0},\phi\bigr{>}_{\Omega}\,\mathrm{d}t={}$$

$$\displaystyle\begin{split}\int_{I}\bigl{<}\partial_{t}v^{0},\psi\bigr{>}_{\Omega}+\langle c^{2}\nabla u^{0},\nabla\psi\rangle_{\Omega}\,\mathrm{d}t={}&\int_{I}\Big{(}\bigl{<}f,\psi\bigr{>}_{\Omega}+\bigl{<}\partial_{n}u,\psi\bigr{>}_{\partial\Omega_{N}}\\[4.30554pt] &\quad-\bigl{<}\partial_{t}v^{D},\psi\bigr{>}_{\Omega}-\langle c^{2}\nabla u^{D},\nabla\psi\rangle_{\Omega}\Big{)}\,\mathrm{d}t\,.\end{split}$$

$$\mathbb{P}_{k}(I_{n};B)=\bigg{\{}w_{\tau}:I_{n}\mapsto B\;\Big{|}\;w_{\tau}(t)=\sum_{j=0}^{k}W^{j}t^{j}\,,\;\forall t\in I_{n}\,,\;W^{j}\in B\;\forall j\bigg{\}}\,.$$

$$X_{\tau }^k(B):=\left\{w_{\tau }\in C(\overline{I};B)\mid w_{\tau }{}_{|I_n}\in {\mathbb{P}}_k(I_n;B)\forall I_n\in {\mathcal{M}}_{\tau }\right\},$$

$$Y_{\tau }^l(B):=\left\{w_{\tau }\in L^2(I;B)\mid w_{\tau }{}_{|I_n}\in {\mathbb{P}}_l(I_n;B)\forall I_n\in {\mathcal{M}}_{\tau }\right\}.$$

$$V_h=V_h^{(p)}=\left\{v_h\in C(\overline{\Omega })\mid v_h{}_{|T}\circ T_K\in {\mathbb{Q}}_p\forall K\in {\mathcal{T}}_h\right\}\cap H_{0,D}^1(\Omega ),$$

$$⟨{\mathcal{A}}_{h}w,v_h⟩=⟨\nabla w,\nabla v_h⟩\text{for all}{v}_h\in V_h.$$

$${\partial }_t^{s_0}{w}_{\tau ,h}^0{}_{|}{}_{I_n}(t_{n-1})={\partial }_t^{s_0}{w}_{\tau ,h}^0{}_{|}{}_{I_{n-1}}(t_{n-1}),\text{for}{w}_{\tau ,h}^0\in \left\{u_{\tau ,h}^0,v_{\tau ,h}^0\right\},$$

$${\partial }_t^{s_1}{u}_{\tau ,h}^0{}_{|}{}_{I_n}(t_n)-{\partial }_t^{s_1-1}{v}_{\tau ,h}^0{}_{|}{}_{I_n}(t_n)=0,$$

$$\begin{array}{rl}{\partial }_t^{s_1}{v}_{\tau ,h}^0{}_{|}{}_{I_n}(t_n)+{\mathcal{A}}_h{\partial }_t^{s_1-1}& u_{\tau ,h}^0{}_{|}{}_{I_n}(t_n)={\partial }_t^{s-1}f(t_n)\\ & -{\partial }_t^{s_1}{v}_{\tau ,h}^D{}_{|}{}_{I_n}(t_n)-{\mathcal{A}}_h{\partial }_t^{s_1-1}{u}_{\tau ,h}^D{}_{|}{}_{I_n}(t_n),\end{array}$$

$$\displaystyle\int_{I_{n}}\Big{(}\langle\partial_{t}u^{0}_{\tau,h},\varphi_{\tau,h}\rangle_{\Omega}-\langle v^{0}_{\tau,h},\varphi_{\tau,h}\rangle_{\Omega}\Big{)}\,\mathrm{d}t=0\,,$$

$$\displaystyle\begin{aligned} \int_{I_{n}}\Big{(}\langle\partial_{t}v^{0}_{\tau,h},\psi_{\tau,h}\rangle_{\Omega}&+\langle\mathcal{A}u^{0}_{\tau,h},\psi_{\tau,h}\rangle_{\Omega}\Big{)}\,\mathrm{d}t=\int_{I_{n}}\langle f,\psi_{\tau,h}\rangle_{\Omega}\,\mathrm{d}t\\ &-\int_{I_{n}}\Big{(}\langle\partial_{t}v^{D}_{\tau,h},\psi_{\tau,h}\rangle_{\Omega}+\langle\mathcal{A}_{h}u^{D}_{\tau,h},\psi_{\tau,h}\rangle_{\Omega}\Big{)}\,\mathrm{d}t\,.\end{aligned}$$

$$\Vert e^u(t)\Vert +\Vert e^v(t)\Vert$$

$$\Vert \nabla e^u(t)\Vert$$

$$\Vert e^u(t){\Vert }_{L^2(I;H)}+\Vert e^v(t){\Vert }_{L^2(I;H)}$$

$$\Vert \nabla e^u(t){\Vert }_{L^2(I;H)}$$

$$\mathbf{M}$$

$$\widehat{{\xi }_0}(0)$$

$$\widehat{{\xi }_1}(0)$$

$$\widehat{{\xi }_2}(0)$$

$$\widehat{{\xi }_3}(0)$$

$$\widehat{{\xi }_0}$$

$$w_{\tau ,h}(\mathbf{x},t)=\sum _{l=0}^3\sum _{j=1}^J{w}_{n,l,j}{\varphi }_j(\mathbf{x}){\xi }_l(t),(\mathbf{x},t)\in \Omega \times \overline{I_n}.$$

$$\mathcal{B}=\{{\varphi }_1{\zeta }_0,\dots ,{\varphi }_J{\zeta }_0\}.$$

$$I_{\tau |I_n}g(t)$$

$$\mathbf{M}\left(-{\mathbf{u}}_{n,0}^0+{\mathbf{u}}_{n,2}^0\right)-{\tau }_n\mathbf{M}\left(\frac{1}{2}{\mathbf{v}}_{n,0}^0+\frac{1}{12}{\mathbf{v}}_{n,1}^0+\frac{1}{2}{\mathbf{v}}_{n,2}^0-\frac{1}{12}{\mathbf{v}}_{n,3}^0\right)=0,$$

$$\begin{array}{rl}& \mathbf{M}\left(-{\mathbf{v}}_{n,0}^0+{\mathbf{v}}_{n,2}^0\right)+{\tau }_n\mathbf{A}\left(\frac{1}{2}{\mathbf{u}}_{n,0}^0+\frac{1}{12}{\mathbf{u}}_{n,1}^0+\frac{1}{2}{\mathbf{u}}_{n,2}^0-\frac{1}{12}{\mathbf{u}}_{n,3}^0\right)=\\ & {\tau }_n\mathbf{M}\left(\frac{1}{2}{\mathbf{f}}_{n,0}+\frac{1}{12}{\mathbf{f}}_{n,1}+\frac{1}{2}{\mathbf{f}}_{n,2}-\frac{1}{12}{\mathbf{f}}_{n,3}\right)-\mathbf{M}\left(-{\mathbf{v}}_{n,0}^D+{\mathbf{v}}_{n,2}^D\right)\\ & -{\tau }_n\mathbf{A}\left(\frac{1}{2}{\mathbf{u}}_{n,0}^D+\frac{1}{12}{\mathbf{u}}_{n,1}^D+\frac{1}{2}{\mathbf{u}}_2^D-\frac{1}{12}{\mathbf{u}}_{n,3}^D\right).\end{array}$$

$${\mathbf{u}}_{n,0}^0={\mathbf{u}}_{n-1,2}^0,{\mathbf{u}}_{n,1}^0={\mathbf{u}}_{n-1,3}^0,{\mathbf{v}}_{n,0}^0={\mathbf{v}}_{n-1,2}^0,{\mathbf{v}}_{n,1}^0={\mathbf{v}}_{n-1,3}^0.$$

$$\mathbf{M}\frac{1}{{\tau }_n}{\mathbf{u}}_{n,3}^0-\mathbf{M}{\mathbf{v}}_{n,2}^0$$

$$\mathbf{M}\frac{1}{{\tau }_n}{\mathbf{v}}_{n,3}^0+\mathbf{A}{\mathbf{u}}_{n,2}^0$$

$$\mathbf{S}\mathbf{x}=\mathbf{b}$$

$$\mathbf{x}={\left({\left({\mathbf{v}}_{n,2}^0\right)}^{\top },{\left({\mathbf{v}}_{n,3}^0\right)}^{\top },{\left({\mathbf{u}}_{n,2}^0\right)}^{\top },{\left({\mathbf{u}}_{n,3}^0\right)}^{\top }\right)}^{\top }$$

Full text

1 Introduction

The accurate and efficient numerical simulation of wave phenomena continues to remain a challenging task and attract researchers’ interest. Wave phenomena are studied in various branches of natural sciences and technology. For instance, fluid-structure interaction, acoustics, poroelasticity, seismics, electro-magnetics and non-destructive material inspection represent prominent fields in that wave propagation is studied. One of our key application for wave propagation is structural health monitoring of lightweight material (for instance, carbon-fibre reinforced polymers) by ultrasonic waves in aerospace engineering. The conceptional idea of this new and intelligent approach is sketched in fig. 1.

The structure is equipped with an integrated actuator-sensor network. The ultrasonic waves that are emitted by the actuators interact with material defects of the solid structure. By means of an inverse modelling, the signals that are recorded by the sensors monitor material failure (cf. [31]) and, as perspective for the future, may allow prognoses about the structure’s residual lifetime. The design of such monitoring systems and the signal interpretation require the elucidation of wave propagation in composite material which demands for highly advanced and efficient numerical simulation techniques.

High-order numerical approximation of partial differential equations has been strongly focused and investigated in the last decades. High order methods are known to be efficient if they approximate functions with large elements of high polynomial degree in regions of high regularity. Prominent examples are $hp$- and spectral element methods in application areas such as computational fluid dynamics or computational mechanics. Their theoretical convergence analysis and adaptive $hp$- and spectral element versions still experience strong development. Whereas high-order approaches have been considered for the approximation of the spatial variables, first- or second-order implicit schemes are often still used for the discretization of the time variable. We note that in this work only implicit time discretization schemes are in the scope of interest. Thus, all remarks refer to this class of methods. Our motivation for using implicit time discretization schemes comes from the overall goal to apply the proposed Galerkin–collocation techniques to mixed systems like, for instance, fluid-structure interaction for free flow modelled by the Navier–Stokes equations [42] or fully dynamic poroelasticity [38].

Driven by the tremendous increase in computing power of modern high performance computing systems and recent progress in the technology of algebraic solver, including efficient techniques of preconditioning, space-time finite element approaches have recently attracted high attention and have been brought to application maturity; cf., e.g., [27, 16, 18, 44]. Space-time finite element methods offer appreciable advantages over discretizations of mixed type based on finite difference techniques for the discretization of the time variable (e.g., by Runge-Kutta methods) and, for instance, finite element methods for the discretization of the space variables. In particular, advantages are the natural embedding of higher order members in the various families of schemes, the applicability of functional analysis techniques in their analyses due to the uniform space-time framework and the applicability of well-known adaptive mesh refinement techniques, including goal-oriented error control [7]. In the meantime a broad variety of implementations of space-time finite element methods does exist. The families of schemes differ by the choices of the trial and test spaces. This leads to continuous or discontinuous approximations of the time variable (cf. e.g., [37, 52]). Further, the fully coupled treatment of all time steps versus time-marching approaches is discussed. In particular, the simultaneous computation of all time steps imposes high demands on the linear solver technology (cf., e.g., [16, 17, 18]).

In this work, we propose the Galerkin–collocation method for the numerical solution of wave equations. This approach combines variational approximation in time by finite element techniques with the concepts of collocation methods and follows the ideas of [13]. By imposing collocation conditions, the test space of the variational condition is downsized. The key ingredients and innovations of the approach are:

• A.

Higher order regularity in time of the fully discrete approximation;

• B.

Linear systems of reduced complexity;

Ingredient [A] is a direct consequence of the construction of the schemes. Higher order regularity is ensured by imposing collocation conditions at the discrete time nodes and endpoints of the subintervals $[t_{n-1},t_{n}]$, for $n=1,\ldots,N$, of the global time interval $[0,T]$. Higher order regularity in time might offer appreciable advantages for future approximations of coupled multi-physics systems if higher order time derivatives of the discrete solution of one subproblem arise as coefficient functions in other subproblems. Ingredient [B] is ensured by the proper choice of a special for the discrete in time function spaces. Thereby, simple vector identities for the degrees of freedom in time are obtained at the left endpoints of the subintervals without generating computational costs. These vector identities can then be exploited to eliminate conditions from the algebraic systems and reduce its size compared to the standard continuous Galerkin–Petrov approximation in time; cf. [12]. In a further work of the authors [4], it is shown that the optimal order of convergence in time (and space) of the underlying finite element discretization is preserved by the Galerkin–collocation approach. The numerical example of section 4.4 that mimics typical studies of structural health monitoring (cf. fig. 1), demonstrates the superiority of the Galerkin–collocation approach over a standard continuous Galerkin–Petrov method admitting continuity and no differentiability in time of the discrete solution.

For the sake of brevity, standard conforming finite element methods are used for the discretization of the spatial variables in this work. This is done since we focus here on time discretization. In the literature it has been mentioned that discontinuous finite element methods in space offer appreciable advantages over continuous ones for the discretization of wave equations; cf., e.g., [10, 20]. The application of, for instance, the symmetric interior penalty disccontinuous Galerkin method (cf. [11, 31] along with a Galerkin–collocation discretization in time, is straightforward.

This work is organized as follows. In section 2 we introduce our prototype model. In section 3 we present its discretization by two families of Galerkin–collocation methods. In section 4, the discrete form of a member of theses families with $C^{1}$-regularity in time is derived. The resulting algebraic system is built and our algebraic solver is described. In section 5 the discrete form of a member of the Galerkin–collocation family with $C^{2}$-regularity in time is derived. For both methods, the results of our numerical experiments are presented and evaluated.

2 Mathematical problem and notation

As a prototype model, we study the wave problem

[TABLE]

In our application of structural health monitoring (cf. fig. 1), $u$ denotes the scalar valued displacement field, $c\in R$ with $c>0$, is a material parameter and $f$ an external force acting on the domain $\Omega\subset\mathbb{R}^{d}$, with $d=2,3$. Further, $g^{u}$ is a prescribed trace on the Dirichlet part $\partial\Omega_{D}$ of the boundary $\partial\Omega=\partial\Omega_{D}\cup\partial\Omega_{N}$, with $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$. By $\partial_{n}$ we denote the normal derivative with outer unit normal vector $\boldsymbol{n}$. Homogeneous Neumann boundary conditions on $\partial\Omega_{N}$ are prescribed for brevity. Finally, $I=(0,T]$ denotes the time domain. Problem (1) is well-posed and admits a unique solution $(u,\partial_{t}u)\in L^{2}(0,T;H^{1}(\Omega))\times L^{2}(0,T;L^{2}(\Omega))$ under appropriate assumptions about the data; cf. [35]. By imbedding, $u\in C([0,T];H^{1}(\Omega))$ and $v\in C([0,T];L^{2}(\Omega))$ is ensured; cf. [36]. Throughout, we tacitly assume that the solution admits all the (improved) regularity being necessary in the arguments.

Our notation is standard. By $H^{m}(\Omega)$ we denote the Sobolev space of $L^{2}(\Omega)$ functions with derivatives up to order $m$ in $L^{2}(\Omega)$. For brevity, we let $H:=L^{2}(\Omega)$ and $V=H^{1}_{0,D}(\Omega)$ be the space of all $H^{1}$-functions with vanishing trace on the Dirichlet part $\partial\Omega_{D}$ of $\partial\Omega$. By $\langle\cdot,\cdot\rangle_{\Omega}$ we denote the inner product in $L^{2}(\Omega)$. For the norms we use $\|\cdot\|:=\|\cdot\|_{L^{2}(\Omega)}$ and $\|\cdot\|_{m}:=\|\cdot\|_{H^{m}(\Omega)}$ for $m\in\mathbb{N}$ and $m\geq 1$. Finally, the expression $a\lesssim b$ stands for the inequality $a\leq C\,b$ with a generic constant $C$ that is indepedent of the size of the space and time meshes.

By $L^{2}(0,T;B)$, $C([0,T];B)$ and $C^{q}([0,T];B)$, for $q\in\mathbb{N}$, we denote the standard Bochner spaces of $B$-valued functions for a Banach space $B$, equipped with their natural norms. Further, for a subinterval $J\subseteq[0,T]$, we will use the notations $L^{2}(J;B)$, $C^{m}(J;B)$ and $C^{0}(J;B):=C(J;B)$ for the corresponding Bochner spaces.

To derive our Galerkin–collocation approach, we first rewrite problem (1) as a first order system in time for the unknowns $(u,v)$, with $v=\partial_{t}u$,

[TABLE]

Further, we represent the unknowns $u$ and $v$ in terms of

[TABLE]

Here, $u^{D}$, $v^{D}\in C(\overline{I};H^{1}(\Omega))$ are supposed to be (extended) functions with traces $u^{D}=g^{u}$ and $v^{D}=g^{v}:=\partial_{t}g^{u}$ on the Dirichlet part $\partial\Omega_{D}$ of $\partial\Omega$.

Using (2) and (3), we then consider solving the following variational problem: Find $(u^{0},v^{0})\in L^{2}(0,T;H^{1}_{0,D}(\Omega))\times L^{2}(0;T;H^{1}_{0,D}(\Omega))$ such that

[TABLE]

and, for all $(\phi,\psi)\in(L^{2}(0;T;H^{1}_{0,D}(\Omega)))^{2}$,

[TABLE]

Remark 1

*i) We note that the correct treatment of inhomogeneous time-dependent boundary conditions is an import issue in the application of variational space-time methods. The space-time discretization that is derived below (cf. section 3) and based on the variational problem (5), (6) ensures convergence rates of optimal order in space and time, also for time-dependent boundary conditions. This is confirmed by the second of the numerical experiments given in section 4.3.

ii) Our Galerkin–collocation approach is based on solving, along with some collocation conditions, the variational equations (5), (6) in finite dimensional subspaces. In particular, the same approximation space will be used for $u^{0}$ and $v^{0}$. For this reason, the solution space for $v^{0}$ and the test space in eq. 6 are chosen slightly stronger than usually; cf. [7]. Choosing $L^{2}(0;T;L^{2}(\Omega)))$ instead, would have been sufficient.*

3 Galerkin-collocation schemes

In this section we introduce two families of Galerkin–collocation schemes. These families combine the concept of collation condition methods applied to the spatially discrete counterpart of the equations (2) with the finite element discretization of the variational equations (5), (6). The collocation constrains then allow us to reduce the size of the discrete test spaces for the variational conditions compared to a standard Galerkin–Petrov approach; cf. [11].

First, we need some notation. For the time discretization we decompose the time interval $I=(0,T]$ into $N$ subintervals $I_{n}=(t_{n-1},t_{n}]$, where $n\in\{1,\ldots,N\}$ and $0=t_{0}<t_{1}<\cdots<t_{n-1}<t_{n}=T$ such that $I=\bigcup_{n=1}^{N}I_{n}$. We put $\tau=\max_{n=1,\ldots N}\tau_{n}$ with $\tau_{n}=t_{n}-t_{n-1}$. Further, the set of time intervals $\mathcal{M}_{\tau}:=\{I_{1},\ldots,I_{n}\}$ is called the time mesh. For a Banach space $B$ and any $k\in\mathbb{N}$, we let

[TABLE]

For an integer $k\in\mathbb{N}$ we introduce the space of globally continuous functions in time

[TABLE]

and for an integer $l\in\mathbb{N}_{0}$ the space of globally $L^{2}$-functions in time

[TABLE]

For the space discretization, let $\mathcal{T}_{h}=\{K\}$ be a shape-regular mesh of $\Omega$ consisting of quadrilateral or hexahedral elements with mesh size $h>0$. Further, for some integer $p\in\mathbb{N}$ let $V_{h}=V_{h}^{(p)}$ be the finite element space that is given by

[TABLE]

where $T_{K}$ is the invertible mapping from the reference cell $\hat{K}$ to the cell $K$ of $\mathcal{T}_{h}$ and $\mathbb{Q}_{p}$ is the space of all polynomials of maximum degree $p$ in each variable. We let $\mathcal{A}_{h}:H^{1}_{0,D}(\Omega)\mapsto V_{h}$ be the discrete operator that is defined by

[TABLE]

Moreover, $(u_{0,h},v_{0,h})\in V_{h}^{2}$ and $(u^{D}_{\tau,h},v^{D}_{\tau,h})\in(C([0,T];V_{h}))^{2}$ define suitable finite element approximations of the initial values $(u_{0},v_{0})$ and the extended boundary values $(u^{D},v^{D})$ in eq. 3. Here, we use interpolation in $V_{h}$ of the given data.

Now we define our classes of Galerkin–collocation schemes. We follow the lines of [4, 13]. We restrict ourselves to the schemes studied in the numerical experiments presented in Secs. 4.3, 4.4 and 5.2. The definition of classes of Galerkin–collocation schemes with even higher regularity in time is straightforward, but not done here.

Definition 1** ($\boldsymbol{C^{l}}$–regular in time Galerkin-collocation schemes GCC$\boldsymbol{{}^{l}(k)}$)**

Let $l\in\{1,2\}$ be fixed and $k\geq 2l+1$. For $n=1,\ldots,N$ and given $(u_{\tau,h}{}_{|}{}_{I_{n-1}}(t_{n-1})$, $v_{\tau,h}{}_{|}{}_{I_{n-1}}(t_{n-1}))\in V_{h}^{2}$ for $n>1$ and $u_{\tau,h}{}_{|}{}_{I_{0}}(t_{0})=u_{0,h}$, $v_{\tau,h}{}_{|}{}_{I_{0}}(t_{0})=v_{0,h}$ for $n=1$, find $(u^{0}_{\tau,h}{}_{|}{}_{I_{n}},v^{0}_{\tau,h}{}_{|}{}_{I_{n}})\in(\mathbb{P}_{k}(I_{n};V_{h}))^{2}$ such that, for $s_{0}\in\mathbb{N}_{0}$, $s_{1}\in\mathbb{N}$ with $s_{0},s_{1}\leq l$,

[TABLE]

and, for all $(\varphi_{\tau,h},\psi_{\tau,h})\in(\mathbb{P}_{0}(I_{n};V_{h}))^{2}$,

[TABLE]

Remark 2

• •

In eq. 11, the discrete initial values $(\partial_{t}u_{\tau,h}(0),\partial_{t}v_{\tau,h}(0))$ arise for $s_{0}=1$. For $\partial_{t}u_{\tau,h}(0)$ we use a suitable finite element approximation $v_{0,h}\in V_{h}$ (here, an interpolation) of $v_{0}\in V$. For $\partial_{t}v_{\tau,h}(0)$ we evaluate the wave equation in the initial time point and use a suitable finite element approximation (here, an interpolation) of $\partial_{t}^{2}u(0)=c^{2}\Delta u(0)+f(0)$. For $s_{0}=2$ in eq. 11, the initial value $\partial_{t}^{2}v_{\tau,h}(0)$ is computed as a suitable finite element approximation (here, an interpolation) of $\partial_{t}^{3}u(0)=c^{2}\Delta\partial_{t}u(0)+\partial_{t}f(0)$. Mathematically, this approach requires that the partial equation and its time derivative are satisfied up to the initial time point and, thereby, sufficient regularity of the continuous solution. Without such regularity assumptions, the application of higher order discretization schemes cannot be justified rigorously. Nevertheless, in practice such methods often show a superiority over lower-order ones, even for solutions without the expected high regularity (cf. section 4.4).

• •

From eq. 11, $(u_{\tau,h},v_{\tau,h})\in(C^{l}(\overline{I};V_{h}))^{2}$, for fixed $l\in\{1,2\}$, is easily concluded.

An optimal order error analysis for the GCC${}^{1}(k)$ family of schemes of Def. 1 is provided in [4]. The following theorem is proved.

Theorem 1** (Error estimates for {\boldsymbol{(u_{\tau,h},v_{\tau,h})}} of

GCC${}^{1}(\boldsymbol{k})$)**

Let $l=1$ and $k\geq 3$. For the error $(e^{\;u},e^{\;v})=(u-u_{\tau,h},v-v_{\tau,h})$ of the fully discrete scheme GCC${}^{l}(k)$ of Def. 1 there holds that

[TABLE]

as well as

[TABLE]

Error estimates for the GCC${}^{2}(k)$ family remain as a work for the future. In section 5.2, the convergence of GCC${}^{2}(5)$ is demonstrated numerically. Further, we note that a computationally cheap post-processing of improved regularity and accuracy for continuous Galerkin–Petrov methods is presented and studied in [12].

In the next sections we study the schemes GCC${}^{1}(3)$ and GCC${}^{2}(5)$ of Def. 1 in detail. Their algebraic forms are derived and the algebraic linear solver are presented. Finally, the results of our numerical experiments with the proposed methods are presented. Here, we restrict ourselves to the lowest-order cases with $k=3$ for $l=1$ and $k=5$ for $l=2$ of Def. 1. This is sufficient to demonstrate the potential of the Galerkin–collocation approach and its superiority over the standard continuous Galerkin approach in space and time [30, 12]. An implementation of GCC${}^{l}(k)$ for higher values of $k$ along with efficient algebraic solvers is currenty still missing.

4 Galerkin–collocation GCC$\boldsymbol{{}^{1}(3)}$

Here, we derive the algebraic system of the GCC${}^{1}(3)$ approach and discuss our algebraic solver for the arising block system. For brevity, the derivation is done for $k=3$ only. The generalization to larger values of $k$ is straightforward; cf. [4].

4.1 Fully discrete system

To derive the discrete counterparts of the variational conditions (14), (15) and the collocation constrains (11) to (13), we let $\{\phi_{j}\}_{j=1}^{J}\subset V_{h}$, denote a (global) nodal Lagrangian basis of $V_{h}$. The mass matrix $\boldsymbol{M}$ and the stiffness matrix $\boldsymbol{A}$ are defind by

[TABLE]

On the reference time interval $\hat{I}=[0,1]$ we define a Hermite-type basis $\{\hat{\xi}_{l}\}_{l=0}^{3}\subset\mathbb{P}_{3}(\hat{I};\mathbb{R})$ of $\mathbb{P}_{3}(\hat{I};\mathbb{R})$ by the conditions

[TABLE]

These conditions then define the basis of $\mathbb{P}_{3}(\hat{I};\mathbb{R})$ by

[TABLE]

By means of the affine transformation $T_{n}(\hat{t}):=t_{n-1}+\tau_{n}\cdot\hat{t}$, with $\hat{t}\in\hat{I}$, from the reference interval $\hat{I}$ to $I_{n}$ such that $t_{n-1}=T_{n}(0)$ and $t_{n}=T_{n}(1)$, the basis $\{\xi_{l}\}_{l=0}^{3}\subset\mathbb{P}_{3}(I_{n};\mathbb{R})$ is given by $\xi_{l}=\hat{\xi}_{l}\circ T_{n}^{-1}$ for $l=0,\ldots,3$. In terms of basis functions, $w_{\tau,h}\in\mathbb{P}_{3}(I_{n};V_{h})$ is thus represented by

[TABLE]

For $\zeta_{0}\equiv 1$ on $\overline{I_{n}}$, a test basis of $\mathbb{P}_{0}(I_{n};V_{h})$ is then given by

[TABLE]

To evaluate the time integrals on the right-hand side of eq. 15 we still use the Hermite-type interpolation operator $I_{\tau|I_{n}}$, on $I_{n}$, defined by

[TABLE]

Here, the values $\partial_{t}^{s}g_{|}{}_{I_{n}}(t_{n-1})$ and $\partial_{t}^{s}g_{|}{}_{I_{n}}(t_{n})$ in (23) denote the corresponding one-sided limits of values $\partial_{t}^{s}g(t)$ from the interior of $I_{n}$.

Now, we can put the equations of the proposed GCC${}^{1}(3)$ approach in their algebraic forms. In the variational equations (14) and (15), we use the representation (21) for each component of $(u_{\tau,h},v_{\tau,h})\in(\mathbb{P}_{3}(I_{n};V_{h}))^{2}$, choose the test functions (22) and interpolate the right-hand side of (15) by applying (23). All of the arising time integrals are evaluated analytically. Then, we can recover the variational conditions (14) and (15) on the subinterval $I_{n}$ in their algebraic forms

[TABLE]

This gives us the first two equations for the set of eight unknown solution vectors $\mathcal{L}=\{\boldsymbol{u}_{n,0}^{0},\ldots,\boldsymbol{u}_{n,3}^{0},\boldsymbol{v}_{n,0}^{0},\ldots,\boldsymbol{v}_{n,3}^{0}\}$ on each subinterval $I_{n}$, where each of these vectors is defined by means of (21) through $\boldsymbol{w}=(w_{1},\ldots,w_{J})^{\top}$ for $\boldsymbol{w}\in\mathcal{L}$.

Next, we study the algebraic forms of the collocations conditions (11) to (13). By means of the definition (19) of the basis of $\mathbb{P}(I_{n};\mathbb{R})$, the constraints (11) read as

[TABLE]

By means of (19) along with (18), the conditions (12) and (13) can be recovered as

[TABLE]

Putting relations (26) into the identities (24) and (LABEL:eq:time_scheme_a2) and combining the resulting equations with (27) and (28) yields for the subinterval $I_{n}$ the linear block system

[TABLE]

for the vector of unknowns

[TABLE]

and the system $\boldsymbol{S}$ and right-hand side $\boldsymbol{b}$ given by

[TABLE]

with \boldsymbol{b}_{n,4}=\boldsymbol{M}\bigl{(}\boldsymbol{v}^{0}_{n,0}+\boldsymbol{v}^{D}_{n,0}-\boldsymbol{v}^{D}_{n,2}+\frac{\tau_{n}}{2}(\boldsymbol{f}_{n,0}+\boldsymbol{f}_{n,2})+\frac{\tau_{n}}{12}(\boldsymbol{f}_{n,1}-\boldsymbol{f}_{n,3})\bigr{)}-\boldsymbol{A}\bigl{(}\frac{\tau_{n}}{2}(\boldsymbol{u}^{0}_{n,0}+\boldsymbol{u}^{D}_{n,0}+\boldsymbol{u}^{D}_{n,2})+\frac{\tau_{n}}{12}(\boldsymbol{u}^{0}_{n,1}+\boldsymbol{u}^{D}_{n,1}-\boldsymbol{u}^{D}_{n,3})\bigr{)}. By means of the collocation constraints (26), the number of unknown coefficient vectors for the discrete solution $(u_{\tau,h}{}_{|I_{n}},v_{\tau,h}{}_{|I_{n}})\in(\mathbb{P}_{3}(I_{n};V_{h}))^{2}$ is thus effectively reduced from eight to four vectors, assembled now in $\boldsymbol{x}$ by (30).

We note that the first two rows of eq. 31 represent the collocation conditions (27) and (28). They have a sparser structure then the last two rows representing the variational conditions which can be advantageous or exploited for the construction of efficient iterative solvers for (29). Compared with a pure variational approach (cf. [24, 25, 30]), more degrees or freedom are obtained directly by computaionally cheap vector identities (cf. (24)) in GCC${}^{l}(k)$ such that they can be eliminated from the overall linear system and, thereby, used to reduce the systems size.

4.2 Solver technology

In the sequel, we present two different iterative approaches for solving the linear system (29) with the non-symmetric matrix $\boldsymbol{S}$. In section 4.4, a runtime comparison between the two concepts is provided. As basic toolbox we use the deal.II finite element library [8] along with the Trilinos library [22] for parallel computations.

4.2.1 1. Approach: Condensing the linear system

The first method for solving (29) is based on the concepts developed in [31]. The key idea is to use Gaussian block elimination within the system matrix $\boldsymbol{S}$ and end up with a linear system with matrix $\boldsymbol{S}_{r}$ of reduced size for one of the subvectors in $\boldsymbol{x}$ in (30) only, and to compute the remaining subvectors of (30) by computationally cheap post-processing steps afterwards. The reduced system matrix $\boldsymbol{S}_{r}$ should have sufficient potenial that an efficient preconditioner for the iterative solution of the reduced system can be constructed. Of course, the Gauss elimination on the block level can be done in different ways. The goal of our approach is to avoid the inversion of the stiffness matrix $\boldsymbol{A}$ in (31) in the computation of the condensed system matrix $\boldsymbol{S}_{r}$ such that a matrix-vector multiplication with $\boldsymbol{S}_{r}$ just involves calculating $\boldsymbol{M}^{-1}$. At least for discontinuous Galerkin methods in space, where $\boldsymbol{M}$ is block diagonal, this is computationally cheap; cf. [11, 31]. We note that a continuous Galerkin approach in space is used here only in order to simplify the notation and since the discretization in time by the combined Galerkin–collocation approach is in the scope of interest.

Here, we choose the subvector $\boldsymbol{u}_{n,2}^{0}$ of $\boldsymbol{x}$ in (30) as the essential unknown, i.e. as the unknown solution vector of the condensed system with matrix $\boldsymbol{S}_{r}$. By block Gaussian elemination we then end up with solving the linear system,

[TABLE]

with right-hand side vector

[TABLE]

The product of $\boldsymbol{A}\boldsymbol{M}^{-1}\boldsymbol{A}$ in (32) mimics the discretization of a fourth order operator due to the appearance of the product of $\boldsymbol{A}$ with its ”weighted” form $\boldsymbol{M}^{-1}\boldsymbol{A}$. Thereby, the conditioning number of the condensed system is strongly increased (cf. [31]) which is the main drawback in this concept of condensing the overall system (29) to (32) for the essential unknown $\boldsymbol{u}_{n,2}^{0}$. On the other hand, since $\boldsymbol{M}$ and $\boldsymbol{A}$ are symmetric and, thus, $\boldsymbol{A}\boldsymbol{M}^{-1}\boldsymbol{A}=(\boldsymbol{A}\boldsymbol{M}^{-1}\boldsymbol{A})^{\top}$, the condensed matrix $\boldsymbol{S}_{r}$ is symmetric such that the preconditioned conjugate gradient method can be applied. Solving systems of type (32) is carefully studied in [11, 31] and the references given therein.

We solve (32) by the conjugate gradient method. The left preconditioning operator

[TABLE]

with positive $\mu\in\mathbb{R}$ chosen such that the spectral norm of $\boldsymbol{P}^{-1}\boldsymbol{S}_{r}$ is minimised, is applied. For details of the choice of the parameter $\mu$, we refer to [11, 31]. Here, we use $\mu=\sqrt{11/2}$. In order to apply the preconditioning operator $\boldsymbol{P}$ in the conjugate gradient iterations, without assembling $\boldsymbol{P}$ explicitly, i.e. to solve the auxiliary system with matrix $\boldsymbol{P}$, we have to solve linear systems for the mass matrix $\boldsymbol{M}$ and the stiffness matrix $\boldsymbol{A}$. For this, we use embedded conjugate gradient iterations combined with an algebraic multigrid preconditioner of the Trilinos library [22]. The overvall algorithm for solving (32) is sketched in fig. 2. The advantage of this approach is that we just have to store $\boldsymbol{M}$ and $\boldsymbol{A}$ as sparse matrices in the computer memory. We never have to assemble the full matrix $\boldsymbol{S}$ from (31), nor do we have to store the reduced matrix $\boldsymbol{S}_{r}$ from (32). Finally, the remaining unknown subvectors $\boldsymbol{v}_{n,2}^{0},\boldsymbol{v}_{n,3}^{0}$ and $\boldsymbol{u}_{n,3}^{0}$ in (29) are successively computed in post-processing steps.

4.2.2 2. Approach: Solving the non-symmetric system

The second approach used to solve (29) relies on assembling the system matrix $\boldsymbol{S}$ of (31) as a sparse matrix and solving the resulting non-symmetric system. For smaller dimensions of $\boldsymbol{S}$ a parallel direct solver [14] is used. For constant time step sizes $\tau_{n}$ the matrix $\boldsymbol{S}$ needs to be factorized once only, which results in excellent performance properties for large sequences of time steps. For high-dimensional problems with interest in practice, we use the Generalized Minimal Residual (GMRES) method, an iterative Krylov subspace method, to solve (29). The drawback of this approach then comes through the necessity to provide an efficient preconditioner, i.e.  an approximation to the inverse of $\boldsymbol{S}$, for the complex block matrix $\boldsymbol{S}$ of (31). Here, we use the algebraic multigrid method as preconditioning technique. We use the MueLue preconditioner [39], which is part of the Trilinos project, with non-symmetric smoothed aggregation. We use an usual V-cycle algorithm along with a symmetric successive over-relaxation (SSOR) smoother with a damping factor of $1.33$. The design of efficient algebraic solvers for block systems like (29), and for higher order variational time discretizations in general, is still an active field of research. We expect further improvement in the future.

4.3 Numerical convergence tests

In this section we present a numerical convergence test for the proposed GCC${}^{1}(3)$ approach of Def. 1 and section 4.1, respectively. For the solution $\{u,v\}$ of eq. 2 and the fully discrete approximation GCC${}^{1}(3)$ of Def. 1 we let

[TABLE]

We study the error $(e^{u},e^{v})$ with respect to the norms

[TABLE]

where $w\in(u,v)$, and in the energy quantities

[TABLE]

All $L^{\infty}$-norms in time are computed on the discrete time grid

[TABLE]

For our first convergence test we prescribe the solution

[TABLE]

on $\Omega\times I=(0,1)^{2}\times[0,1]$. We let $c=1$, use a constant mesh size $h_{0}=0.25$ and start with the time step size $\tau_{0}=0.1$. We compute the errors on a sequence of successively refined time meshes by halving the step sizes in each refinement step. We choose a bicubic discretization of the space variables in $V_{h}^{3}$ (cf. (9)) such that the spatial part of the solution is resolved exactly by its numerical approximation. Table 1 summarizes the computed errors and experimental orders of convergence. The expected convergence rates of Thm. 1 are nicely confirmed.

In our second numerical experiment we study the space-time convergence behavior of a solution satisfying non-homogeneous Dirichlet boundary conditions,

[TABLE]

on $\Omega\times I=(0,1)^{2}\times[0,1]$. We choose a bicubic discretization in $V_{h}^{3}$ (cf. (9)) of the space variable. We refine the space-time mesh by halving both step sizes in each refinement step. Table 2 shows the computed errors and experimental orders of convergence for this example. In all measured norms, optimal rates in space and time (cf. Thm. 1) are confirmed. This underlines the correct treatment of the prescribed non-homogeneous Dirichlet boundary conditions.

4.4 Test case of structural health monitoring

Next, we consider a test problem that is based on [7] and related to typical problems of structural health monitoring by ultrasonic waves (cf. fig. 1). We aim to compare the GCC${}^{1}(3)$ approach with a standard continuous in time Galerkin–Petrov approach cGP(2) of piecewise quadratic polynomials in time; cf. [12, 24] for details. The cGP(2) scheme has superconvergence properties in the discrete time nodes $t_{n}$ for $n=1,\ldots,N$ as shown in [12]. Thus,the errors $\max_{n=1,\ldots,N}\|e^{u}(t_{n})\|$ and $\max_{n=1,\ldots,N}\|e^{v}(t_{n})\|$ for the GCC${}^{1}(3)$ and the cGP(2) scheme admit the same fourth order rate of convergence in time and, thus, are comparable with respect to accuracy.

The test setting is sketched in fig. 3(a). We consider $\Omega\times I=(-1,1)^{2}\times(0,1)$, let $f=0$ and, for simplicity, prescribe homogeneous Dirichlet boundary conditions such that $u^{D}=0$. For the initial value we prescribe a regularized Dirac impulse by

[TABLE]

where $\Theta$ is the Heaviside function. The coefficient function $c(\boldsymbol{x})$, mimicing a material parameter, has a jump discontinuity and is given by $c(\boldsymbol{x})=1$ for $x_{2}<0.2$ and $c(\boldsymbol{x})=9$ for $x_{2}\geq 0.2$. Further we put $v_{0}=0$ for the second initial value. Finally, we define the control region $\Omega_{c}=(0.75-h_{c},0.75+h_{c})\times(-h_{c},h_{c})$ where we calculate the signal arrival, at a sensor position for instance, in terms of

[TABLE]

We choose a spatial mesh of $65\,536$ cells and $\mathbb{Q}_{7}$ elements; cf. fig. 3(b). This leads to more than $3.2\text{\times}{10}^{6}$ degrees of freedom in space in each time step for each of the solution vectors. For each computation of the control quantity (42), with $t\in(0,1]$, we use a constant time step size $\tau_{n}$ for all time steps and compare the computation with the initially chosen reference time step size of $\tau_{0}=$2\text{\times}{10}^{-5}$$.

Figure 4 shows the signal arrival and control quantity (42) over $t\in(0.6,1)$ with different choices of the time step sizes for the Galerkin–collocation scheme GCC${}^{1}(3)$ and the standard Galerkin–Petrov approach cGP(2) (cf. [12, 24]) of a continuous in time approximation. For the cGP(2) approach, very small time step sizes are required to avoid over- and undershoots in the control quantity $u_{c}(t)$. For the GCC${}^{1}(3)$ approach with $C^{1}$ regularity in time, much larger time steps, approximately 100 times larger, can be applied without loss of accuracy compared to the fully converged reference solution given by GCC${}^{1}()3$ with step size $\tau_{0}$. This clearly shows the superiority of the Galerkin–collocation scheme GCC${}^{1}(3)$.

In Table 3 the computational costs are summarized, where $r_{1}$ is the runtime for solving the condensed system by the approach of section 4.2.1 and $r_{2}$ is the runtime for solving the block system by the approach of section 4.2.2. For the cGP(2) approach, only the first of the either iterative solver techniques was implemented. Recalling from fig. 4 that GCC${}^{1}(3)$ with $\tau_{n}=100\times\tau_{0}$ leads to the fully converged solution whereas cGP(2) with $\tau_{n}=\tau_{0}$ already shows over- and undershoots, a strong superiority of GCC${}^{1}(3)$ over cGP(2) is observed in in Table 3. For both solver, a reduction in the wall clock time by a factor of about 25 is shown.

5 Galerkin–collocation GCC$\boldsymbol{{}^{2}(5)}$

Here, we briefly derive the algebraic form of the Galerkin–collocation scheme GCC${}^{2}(k)$ of Def. 1 with fully discrete solutions $(u_{\tau,h}{}_{|I_{n}},v_{\tau,h}{}_{|I_{n}})\in(X_{\tau}^{k}(V_{h}))^{2}$ such that $(u_{\tau,h},v_{\tau,h})\in(C^{2}(\overline{I};V_{h}))^{2}$. For brevity, we restrict ourselves to the lowest polynomial degree in time $k=5$ that is possible to get $C^{2}$-regularity. The convergence properties are then demonstrated numerically.

5.1 Fully discrete system

We follow the lines of section 4.1 and use the notation introduced there. The six basis function of $\mathbb{P}_{5}(\hat{I};\mathbb{R})$ on the reference interval $\hat{I}$ are defined by the conditions

[TABLE]

where $\delta_{i,j}$ denotes the usual Kronecker symbol. This gives us

[TABLE]

For this basis of $\mathbb{P}_{5}(\hat{I};\mathbb{R})$, the discrete variational conditions (14), (15) then read as

[TABLE]

In the basis, the first collocation conditions (11) yield for $\boldsymbol{w}_{n,i}^{0}\in\{\boldsymbol{u}_{n,i}^{0},\boldsymbol{v}_{n,i}^{0}\}$ that

[TABLE]

which reduces the number of unknown solution vectors by $6$ on each subinterval $I_{n}$. For the collocation conditions (12), (13) at $t_{n}$ and $s=1$ we deduce that

[TABLE]

Similarly, for $s=2$ the collocation conditions (12), (13) at $t_{n}$ read as

[TABLE]

Finally, we recover the previous conditions as the linear system $\boldsymbol{S}\boldsymbol{x}=\boldsymbol{b}$ for the vector of unknowns $\boldsymbol{x}=\left((\boldsymbol{u}_{n,3}^{0})^{\top},(\boldsymbol{u}_{n,4}^{0})^{\top},(\boldsymbol{v}_{n,5}^{0})^{\top},\boldsymbol{u}_{n,5}^{0})^{\top}\right)^{\top}$ and with the system matrix $\boldsymbol{S}$ and right-hand side vector $\boldsymbol{b}$ given by

[TABLE]

with \boldsymbol{b}_{n,3}=\boldsymbol{M}\bigl{(}\boldsymbol{u}_{n,0}^{0}+\boldsymbol{u}_{n,0}^{D}-\boldsymbol{u}_{n,3}^{D}+\frac{\tau_{n}}{2}(\boldsymbol{v}_{n,0}^{0}+\boldsymbol{v}_{n,0}^{D})+\frac{\tau_{n}}{10}(\boldsymbol{v}_{n,1}^{0}+\boldsymbol{v}_{n,1}^{D})+\frac{\tau_{n}}{120}(\boldsymbol{v}_{n,2}^{0}+\boldsymbol{v}_{n,2}^{D})+\tau_{n}(\frac{1}{2}\boldsymbol{v}_{n,3}^{D}-\frac{1}{10}\boldsymbol{v}_{n,4}^{D}+\frac{1}{120}\boldsymbol{v}_{n,5}^{D}) and \boldsymbol{b}_{n,4}=\boldsymbol{M}\bigl{(}\boldsymbol{v}_{n,0}^{0}+\boldsymbol{v}_{n,0}^{D}-\boldsymbol{v}_{n,3}^{D}\bigr{)}+\tau_{n}\bigl{(}\frac{1}{2}\boldsymbol{f}_{n,3}+\frac{1}{10}\boldsymbol{f}_{n,1}+\frac{1}{120}\boldsymbol{f}_{n,2}+\frac{1}{2}\boldsymbol{f}_{n,3}-\frac{1}{10}\boldsymbol{f}_{n,4}+\frac{1}{120}\boldsymbol{f}_{n,5}\bigr{)}-\tau_{n}\boldsymbol{A}\bigl{(}\frac{1}{2}(\boldsymbol{u}_{n,0}^{0}+\boldsymbol{u}_{n,0}^{D})2+\frac{1}{10}(\boldsymbol{u}_{n,1}^{0}+\boldsymbol{u}_{n,1}^{D})+\frac{1}{120}(\boldsymbol{u}_{n,2}^{0}+\boldsymbol{u}_{n,2}^{D})+\frac{1}{2}\boldsymbol{u}_{n,3}^{D}-\frac{1}{10}\boldsymbol{u}_{n,4}^{D}+\frac{1}{120}\boldsymbol{u}_{n,5}^{D}\bigr{)}.

5.2 Iterative solver and convergence study

To solve the linear system $\boldsymbol{S}\boldsymbol{x}=\boldsymbol{b}$ with $\boldsymbol{S}$ from (53), we use block Gaussian elimination, as sketched in section 4.2.1, to find a reduced system $\boldsymbol{S}_{r}\boldsymbol{u}_{n,4}^{0}=\boldsymbol{b}_{r}$ for the essential unknown $\boldsymbol{u}_{n,4}^{0}$. All remaining unknown subvectors of $\boldsymbol{x}$ can be computed in post-processing steps. In explicit form, the condensed system reads as

[TABLE]

For brevity, we omit the exact definition of $\boldsymbol{b}_{r}$ that can be deduced easily from (53).

The matrix $\boldsymbol{S}_{r}$ is symmetric such that preconditioned conjugate gradient iterations are used for its solution. The preconditioner is constructed along the lines of section 4.2.1. The remainder part $\tau_{n}^{6}\boldsymbol{A}\boldsymbol{M}^{-^{1}}\boldsymbol{A}\boldsymbol{M}^{-^{1}}\boldsymbol{A}$ is still ignored in the construction of the preconditioner. Even though the remainder is weighted by the small factor $\tau_{n}^{6}$, numerical experiments indicate that this scaling is not sufficient to balance its impact on the interation process. For the construction of an efficient preconditioning technique for $\boldsymbol{S}_{r}$ of GCC${}^{2}(5)$ further improvements are still necessary.

To illustrate the convergence behavior and performance of the GCC${}^{2}(5)$ Galerkin–collocation approach, we present in Table 4 our numerical results for the test problem (39). The expected convergence of sixth order in time is nicely observed in all norms.