Numerical methods for changing type systems

In this note we develop a numerical method for partial differential equations with changing type. Our method is based on a unified solution theory found by Rainer Picard for several linear equations from mathematical physics. Parallel to the solution theory already developed, we frame our numerical method in a discontinuous Galerkin approach in space-time with certain exponentially weighted spaces.


Introduction
Following the rationale presented in [9], most of the classical linear partial differential equations arising in mathematical physics share a common form, namely the form of an evolutionary problem. That is, we consider equations of the form where F is a given source term, ∂ t stands for the derivative with respect to time, M 0 , M 1 are bounded linear operators on some Hilbert space H and A is an unbounded skewselfadjoint operator in H. We are seeking for a unique solution U of the above equation.
We remark here that we do not impose initial conditions, since we consider the whole real line as time horizon, and hence, we implicitly assume a vanishing initial value at "−∞". To illustrate the setting, we begin with presenting some examples.
Example 1.1. Let Ω ⊆ R n an open non-empty set, where n ∈ N, but, typically n ∈ {1, 2, 3}. We define the following two differential operators assigning each function u ∈ H 1 0 (Ω) its gradient, that is, the column-vector of its partial derivatives in each coordinate direction. Moreover, we set div := −(∇ 0 ) * : D(div) ⊆ L 2 (Ω) n → L 2 (Ω), which is nothing but the operator assigning each L 2 vector-field its distributional divergence with maximal domain, that is, Since both the operators ∇ 0 and div are closed and skew-adjoints of one another, we infer that the operator is skew-selfadjoint on the Hilbert space H = L 2 (Ω) × L 2 (Ω) n . Choosing M 0 = 1 and M 1 = 0 in (1.1), the corresponding evolutionary problem reads as If g = 0, this is nothing but the wave equation. Indeed, the second line then gives ∂ t v = −∇ 0 u, and hence, differentiating the first line with respect to time, we obtain Note that div ∇ 0 = ∆ D is the classical Dirichlet-Laplace operator on L 2 (Ω).
Choosing M 0 = 1 0 0 0 and M 1 = 0 0 0 1 in (1.1), the corresponding problem reads as Setting again g = 0, the latter gives the heat equation. Indeed, the second line reads v = −∇ 0 u and hence the first line yields

November 9, 2016
Finally, choosing M 0 = 0 and M 1 = 1 in (1.1), we get which in the case g = 0 gives the elliptic equation Remark 1.2. We note that we can treat the case of homogeneous Neumann boundary conditions in the same way. The only difference is that we define ∇ as the distributional gradient on H 1 (Ω) and div 0 := −(∇) * . Replacing now ∇ 0 by ∇ and div by div 0 yields the same hyperbolic, parabolic and elliptic type problem above, but now with homogeneous Neumann boundary conditions. Example 1.1 shows that evolutionary problems cover all three classical types of partial differential equations, elliptic, parabolic and hyperbolic. However, also problems of mixed type are covered as the next example shows. χ Ω ell 0 0 χ Ωpar∪Ω ell . The resulting evolutionary problem then is of mixed type. More precisely, on Ω ell we get an equation of elliptic type, on Ω par the equations becomes parabolic while on Ω hyp the problem is hyperbolic.
Remark 1.4. The interested reader might wonder that there is not imposed any transmission condition on the unknown quantities along the interfaces of Ω ell , Ω par and Ω hyp . However, this can be implemented automatically by being in the domain of the corresponding operator sum, as can be seen, for instance, in [23,Remark 3.2], see also [13,An illustrative Example]. Another example of a mixeed tyoe problem in contral theory can be found in [12,Remark 6.2] In [9], the well-posedness of problems of the form (1.1) has been addressed. In fact, it was shown that these probolems also cover the classical Maxwell's equations, the equations of linearized elasticity or a general class of coupled phenomena, see, for instance, [7,8,11]. All these problems are indeed well-posed (see Section 2 for the precise statement). The purpose of the present article is to provide numerical methods for such problems. In this article, for the applications to follow, we will focus, however, on problems of mixed type of the form sketched in Example 1.3. Moreover, as the spatial discretisation has to be developed for each problem separately, anyway, in this work, we will put an emphasize on the time-discretisation. Furthermore, we want to stress that the null-space of M 0 in (1.1) might be infinite-dimensional. Hence, we seek to develop a numerical scheme, which in particular allows for the treatment of a certain class of (partial) differential-algebraic equations.

November 9, 2016
For the numerical treatment of the time derivatives we use a discontinuous Galerkin (dG) method, see also Section 3. The first dG-method was published in 1973 on neutron transport [15]. Later the methodology was developed further for classical hyperbolic, parabolic and elliptic problems, see also the survey article [4] and the book [16]. Note that there is a strong connection between dG-methods and Runge-Kutta (collocation) methods, see [2] for parabolic problems. In Section 2, for convenience, we will recall some essentials for evolutionary equations. In particular, we recall the solution theory of problems of the type of equation (1.1). We will introduce a semi-dicretised version, Equation (3.1), of equation (1.1) at the beginning of Section 3. We will also provide a solution theory for this semi-discretised variant with general underlying (spatial) Hilbert space H (Proposition 3.1). The remainder of Section 3 is devoted to estimate difference of the exact solution of (1.1) and the approximate solution of (3.1): In Subsection 3.1, we bound the error by solely in terms of the interpolation error, which will eventually be estimated in Subsection 3.2. As our prime example, we address the full space-time discretisation of Example 1.3 and derive corresponding error estimates. We verify our theoretical findings in Section 5 by means of a 1 + 1-and a 1 + 2-dimensional numerical example. This article is attached an appendix (Section 6), where, for the convenience of the reader, we recall some well-known results on the Gauß-Radau quadrature rule including the fact that the choice of Gauß-Radau points depends continuously on the weighting function. We will need some implications of the fact just mentioned in our a-priori analysis in Subsection 3.1.

The setting of evolutionary problems
In this section we briefly recall the well-posedness result stated in [9]. For doing so, we need to specify the functional analytic setting. Throughout, let H be a real Hilbert space.
Definition. Let ρ > 0 and define the space where we as usual identify functions which are equal almost everywhere. The space H ρ (R; H) is a Hilbert space endowed with the natural inner product given by Moreover, we define ∂ t to be the closure of the operator where by C ∞ c (R; H) we denote the space of infinitely differentiable H-valued functions on R with compact support. We denote the domain of ∂ k t by H k ρ (R; H) for k ∈ N.
Within the setting introduced, we can formulate the well-posedness for evolutionary equations of the form (1.1).
Then, for each ρ ≥ ρ 0 and each F ∈ H ρ (R; H) there exists a unique U ∈ H ρ (R; H) such that where the closure is taken in H ρ (R; H). Moreover, the following continuity estimate holds If F ∈ H k ρ (R; H) for k ∈ N, then so is U and we can omit the closure bar in (2.1).

Remark 2.2.
(a) Note that the positive definiteness condition in the latter theorem especially implies which yields that U (t) ∈ D(A) for almost every t ∈ R. If even F, U ∈ H 2 ρ (R; H) the latter gives AU ∈ H 1 ρ (R; H) and hence, using the Sobolev embedding result (see part (b)), U ∈ C ρ (R; D(A)).
(d) The original result in [9] treat a general class of time-translation invariant coefficients. We refer to [13,22] for non-autonomous variants as well as to [19,20] for non-autonomous and/or non-linear versions of Theorem 2.1.
We note that the equations treated in Example 1.1 and Example 1.3 satisfy the conditions of the previous theorem and hence, are well-posed.

Semi-discretisation in time
In this section, we discretise (1.1) with respect to time and do the a-priori analysis. We assume that holds for all polynomials p of degree at most 2q. Note that the weights and nodes can always be numerically computed as shown for instance in [14,Chapter 4.6], see also the appendix (Section 6) for some basic facts on the Gauß-Radau quadrature. With the following standard transformation we define by instead of the scalar products a, b ρ we employ the following discrete quadrature formulation: For given F ∈ U τ and x 0 ∈ H, find U ∈ U τ , such that for all Φ ∈ U τ and m ∈ {1, 2, . . . , M } it holds Here, we denote by and by Φ + m−1 := Φ(t m−1 +).
Then there exists a unique solution of (3.1).
Proof. Let m ∈ {1, . . . , M } and recall that P q (I m ; H) is a Hilbert space with the aforementioned scalar product. We note that Moreover, the mapping Ψ : H → P q (I m ; H) is linear and bounded, since We now prove that for each g ∈ P q (I m ; H) there is a unique u ∈ P q (I m ; D(A)) such that

November 9, 2016
For doing so, we first compute using integration by parts Therefore, for all u ∈ P q (I m ; D(A)) we get where we have used In particular, both B := (∂ t M 0 + M 1 ) + ΨM 0 δ m−1 and B + A are strictly positive definite. Moreover, since B is bounded, B * is strictly positive definite, as well. Hence, from we read off that (B + A) * is strictly positive definite as well. Thus, for each g ∈ P q (I m ; H) there is a unique u ∈ P q (I m ; D(A)) = D(A + B) such that We put U | Im := u. The thus defined function U solves (3.1): We observe by definition for all Φ ∈ U τ and m ∈ {1, . . . , M }. The latter is the same as saying But, since the quadrature is exact for polynomials up to degree 2q, the latter equation in turn is equivalent to which yields existence of U . Uniqueness follows from the uniqueness of u satisfying (3.2).

On some a-priori error estimates in time
After having proved the unique solvability of (3.1), we address the error estimates in the following. In our analysis we will use the discretised norms as approximations of |||v||| 2 ρ,m := Im |v(t)| 2 H exp(−2ρ(t − t m−1 )) dt and |v| 2 ρ . Note that for v ∈ U τ the approximation is exact. Let us start by defining an interpolation operator into U τ and define by ϕ m,i with i ∈ {0, . . . , q} the associated Lagrange basis functions to the nodes t m,i . Then we obtain for a function v ∈ C([0, T ], H) by an interpolation operator in time.
In the analysis to follow, we will consider the problem (2.1). In particular, we emphasize that we assume that the hypotheses of Theorem 2.1 are in effect.
Furthermore, we fix a right-hand side
We consider the following splitting Note that for almost every t ∈ [0, T ] we have that and thus, since the P F is interpolates at the Gauß-Radau points used in the quadrature. Hence, U solves the same semi-discretised problem as U τ . Thus, we obtain with χ ∈ U τ as test function the error equation for all m ∈ {1, . . . , M }, where the subscripts d and i should remind of discretisation and interpolation, respectively.

November 9, 2016
Proof. Let m ∈ {1, . . . , M }. Since ξ is a (piece-wise) polynomial of order q in time, we obtain Further, we compute Together with the lemma is proved.
In order to analyse E m i we introduce another interpolation operator, that enables us to estimate the time derivative of the interpolation error with a higher order. This operator utilises t m,−1 := t m−1 in addition to t m,i , i ∈ {0, . . . , q} as interpolation points. Denoting the associated Lagrange basis functions by ψ m,i , i ∈ {−1, 0, . . . , q}, this interpolation operator is given by Note the P maps to functions that are continuous in time (recall that t m,q = t m ) while the image of P is allowed to be discontinuous at the time mesh points.
Proof. With U being continuous in time, we only have to consider the discrete part.
Combining the previous lemmas gives a first result.
Theorem 3.5. There exists a C ≥ 0 depending on T , ρ and γ, only, such that Proof. Combining Lemmas 3.2 to 3.4 for ψ = ξ we have for some C ≥ 1 depending on T and ρ only by ξ − 0 = 0 and neglecting the positive jump-contributions, and Thus for β < γ/2 the result is proved upon the equality Remark 3.6. Let m ∈ {1, . . . , M }. Note that the estimate in Theorem 3.5 remains valid, if one respectively replaces T by t m , ξ − M by ξ − m as well as the |||·||| Q,ρ by |||·χĨ||| Q,ρ with χĨ being the characteristic function ofĨ = m k=1 I k .
In the following, we want to improve Theorem 3.5. In order to do so, we will need the following technical lemmas.
and the corresponding integration by parts formula . . , q} be the points and weights of the right-sided Gauß-Radau quadrature rule of order q on (0, 1] with weighting function t → e −2ρt . Then Thus, where we denote by m the multiplication-with-the-argument, that is, (mf )(t) := tf (t). With (3.8) we obtain for the second term From mv Λ ∈ P 2q−1 and mΛ Λ ∈ P 2q together with the exactness of the quadrature rule it follows that and in the same way Λ , mΛ ρ = 2ρ Λ, 1 ρ + e −2ρ − Λ(0). (3.9) We have thus this far Next, (3.9) yields Using Λ, 1 − mΛ ρ ≥ 0, which we provide in Lemma 3.8, the first result is proved. The second one follows upon using the exactness of the quadrature rule and t −1 i > 1.
Proof. We rewrite the scalar product as a quadrature error: for f given by f (t) = tΛ 2 (t), where Q[g] = q i=0 w i g(t i ) and I[g] = 1 0 g for suitable g. There exists a constant α ∈ R and a polynomial w 0 ∈ P q−1 [0, 1], such that where w 1 ∈ P 2q [0, 1]. Thus, setting g(t) = t 2q+1 , we have that due to the exactness of the quadrature rule for polynomials of degree 2q. Let Πw ∈ P 2q [0, 1] be an Hermite-interpolant of a given function w satisfying Then it follows Using that for each t ∈ [0, 1] there is ζ ∈ (0, 1) such that see, for instance, [18, Section 2.1.5], we infer that Now we are able to improve Theorem 3.5 following [1, Corollary 2.1] and [21].
Theorem 3.9. There exists C ≥ 0 depending on T , q, M 0 , M 1 , γ, and ρ such that with g(U ) defined as in Theorem 3.5.
Proof. For the discrete error ξ = U τ − P U ∈ U τ we define ϕ by Then for all m ∈ {1, . . . , M } and i ∈ {0, . . . , q} we have and by Lemma 3.7 (apply the lemma to the functions p = By the equivalence of norms on P q ([0, 1]), there exists K 1 ≥ 0 depending on q only, such that Consequently, we obtain for all m ∈ {1, . . . , M }

Estimating the interpolation error in time
In the previous section we showed that the discrete error is bounded in terms of the interpolation errors. We finalize the error estimates in time in this section focussing on the interpolation error. The aim and, thus, main theorem of this section is Theorem 3.13, where we estimate the difference between the exact solution U of (3.4) and the solution U τ of the quadrature formulation (3.1) with right-hand side P F and initial value U (0+). We use the same notation as in the previous section. In addition, we set τ := max{τ m : m ∈ {1, . . . , M }}. Moreover, shall further assume that the hypotheses of Theorem 2.1 are in effect.
Proof. First we note that H q+3 ρ (R; H) → C q+2 ρ (R; H) by the Sobolev-embedding theorem. By the definition of |||·||| Q,ρ we have that Using the standard result from interpolation theory For the next two lemmas, we recall the standard result from interpolation theory Proof. We obtain The claim follows from the Sobolev-embedding theorem.
With the previous lemmas we can already estimate g(U ). Now let us estimate the final term needed to estimate the error U − U τ .
Proof. Using the Cauchy-Schwarz and Young inequality we derive Using  Combining these results the claim follows from the Sobolev-embedding theorem.
Combining the previous lemmas, Theorem 3.5 and Theorem 3.9, we can bound the discrete error in time.
Remark 3.14. The above analysis holds for all evolutionary problems and gives error bounds for the semi-discrete solution of order q + 1, assuming enough regularity in time.
For a fully discrete method, a spatial discretisation has to be defined too. This step, however, has to be done for each problem considered separately.

Full discretisation for Example 1.3
Let us assume a regular, quasi uniform and shape-regular triangulation Ω h of Ω into triangular open cells σ with maximal cell diameter h. Moreover, we assume that the interfaces between Ω ell , Ω par and Ω hyp are polygonal such that the triangulation Ω h fits to these interfaces. As the whole article is mainly concerned with the correct time-discretization, in this section, we will employ the custom of the "generic constant" C ≥ 0 that may vary from line to line, which, however, depends on T , ρ, M 1 , M 0 , q, and γ and on k, the order of the assumed spatial regularity, only. Then the fully discretised counterpart U τ h to U is given by , v h | Im ∈ P q (I m , V 2 (Ω)), m ∈ {1, . . . , m}} , where the spatial spaces are V 1 (Ω) := v ∈ H 1 0 (Ω); ∀σ : v| σ ∈ P k (σ) , V 2 (Ω) := {w ∈ H(div, Ω); ∀σ : w| σ ∈ RT k−1 (σ)} .
Here P k (σ) is the space of polynomials of degree up to k on the cell σ and RT k−1 (σ) is a the Raviart-Thomas-space, defined by

Note that
Furthermore, if the mesh consists of quadrilateral or hexahedral cells, in above definitions and statements the polynomials space P k (σ) can be replaced by Q k (σ) including all polynomials of total degree k over σ.
Remark 4.1 (Solvability of the fully discrete system). We can apply the general existence theory that was also used in Proposition 3.1. More precisely, the positive definiteness still holds, since the triangulation fits to the interfaces and hence, the uniqueness of the system is warranted. However, since the problem is finite-dimensional, the uniqueness implies the existence of a solution of the fully discretised problem.
Let us come to the interpolation operator I = (I 1 , I 2 ). For I 1 : C(Ω) → V 1 we use the Scott-Zhang interpolant on each cell σ, see [17] for a precise definition, that is patched together continuously. Here local interpolation error estimates can be given using L 2norms also in 3d, which is not possible for standard Lagrange interpolation. For I 2 : we also use the standard interpolator, defined via moments, see [3]. Note that in the following, in order to avoid a cluttered notation as much as possible, we will not explicitly keep track on the number of components of the L 2 (Ω)-or H k (Ω)-spaces under consideration, as it will be obvious from the context.

Standard local interpolation error estimates yield for all
where 1 ≤ r ≤ k + 1 and for all q ∈ H s (Ω) such that div q ∈ H s (Ω) q − I 2 q 0,Ω ≤ Ch s q s,Ω , div(q − I 2 q) 0,Ω ≤ Ch s div q s,Ω , where 1 ≤ s ≤ k, see [3]. Let U τ h ∈ U τ h be the solution of the fully discretised system and P IU ∈ U τ h be the interpolated solution of (1.1) for the operators M 0 , M 1 given in Example 1.3 and A given as in Example 1.1. Then we obtain analogously to the derivation of the errors of the semi-discretisation where we remark that in contrast to Theorem 3.5 the terms |||M 1 (U − P IU )||| 2 Q,ρ and |||A(U − P IU )||| 2 Q,ρ do not vanish, since we also interpolate with respect to space. In the following group of lemmas we estimate the terms on the right-hand side of (4.1) and start with a term partcularly needed for the final convergence estimate in Theorem 4.7. Beforehand, let us introduce where D ⊆ Ω is measurable.
Proof. By the definition of Q m [·] ρ we have Very similarly we have for the second norm Proof. The assertion follows from Lemma 4.2 and the boundedness of M 1 .
Proof. The operator M 0 is selfadjoint and non-negative. Thus it follows that for each t ∈ [0, T ]. The second term can be estimated by according to Lemma 3.12, while the first term can be estimated by due to the boundedness of √ M 0 . Hence, the assertion follows.
Proof. We have that by Lemma 3.10. For the first term we have by Lemma 4.2 Proof. This is a direct consequence of Lemma 3.11. Theorem 4.7. We assume for the solution U = (u, v) of Example 1.3 the regularity as well as AU ∈ H ρ (R; H k (Ω) × H k (Ω)).
Then we have for the error of the numerical solution by (3.1)

Numerical examples
In the following section we consider some examples to verify numerically our theoretical findings.
Again the exact solution can be found and is given by . Note that u and v are non-differentiable, but piece-wise smooth. Figure 3 shows the solutions for t ∈ [0, 1]. Note that a priori, we impose no transmission condition. However, as in [23,Remark 3.2], they can be derived for u satisfying (5. we consider in Table 3 we observe a convergence behaviour similar to the previous smooth case.    , Ω e = Ω \Ω h and Ω p = ∅. The problem is given on (0, T ) × Ω by where f (t, x) = 2 sin(πt)χ R <1/2 ×R (x).   shows some snapshots of the component u of the solution U , approximated by a numerical simulation. In order to investigate the error-behaviour upon refinement of the discretisation, we use a numerically computed reference solutionŨ instead of the real one U . For this we set T = 1 and use an equidistant mesh of 128×128 rectangular cells in space and 128 cells in time, and polynomial degrees p = 3 and q = 2. Thus u is approximated in space by piece-wise Q 3 elements, v by RT 2 -elements and both in time by P 2 -elements. In Table 5 we see the results of our numerical simulation for two pairs of polynomial order. We observe, that the error rates are independent of the polynomial order and furthermore less than the optimal orders given in Theorem 4.7. The reason for this decrease in convergence order lies in the reduced regularity of the solution to this given problem. The interior boundaries where the type of the problem changes introduces corners, where it is very likely for singular solution components to arise.
Then, using that f, p q w = 0, we compute x − r k r j − r k w(x)dx = q j=0 g(r j )ω j .
We address the continuous dependence of the Gauß-Radau points on the weighting function as follows.
Then, for every compact set K ⊂ R, we have Proof. Assume by contradiction that there exists (τ n ) n convergent to some τ with the property But, which contradicts the assumption.
The next two corollaries are the ones needed in Subsection 3.1, that is, for the error estimate with respect to the time-discretization. Beforehand, we introduce for a bounded interval I ⊆ R the mapping ϕ I : (−1, 1) → I, where a := inf I, b := sup I. Further, we set |I| := b − a. Hence, the assertion follows from Corollary 6.4.
The next corollary is concerned with the lowest Gauß-Radau point for different weights: Corollary 6.6. For τ ∈ R let w τ and (ω (τ ) , r (τ ) ) be given as in Corollary 6.4. Let T > 0.
Then there exists c > 0 such that for all intervals I ⊆ R with |I| ≤ T and 0 ≤ τ ≤ T we have ϕ I (r (τ ) 0 ) − inf I ≥ c|I|.