Heiko Gimperlein
email: h.gimperlein@hw.ac.uk.

Ceyhun Özdemir

David Stark

Ernst P Stephan

A residual a posteriori error estimate for the time-domain boundary element method

Keywords: Mathematics Subject Classification: 65N38 (primary), 65M15, 35L67 (secondary) boundary element method, a posteriori error estimates, adaptive mesh refinements, screen problems, wave equation

This article investigates residual a posteriori error estimates and adaptive mesh refinements for time-dependent boundary element methods for the wave equation. We obtain reliable estimates for Dirichlet and acoustic boundary conditions which hold for a large class of discretizations. Efficiency of the error estimate is shown for a natural discretization of low order. Numerical examples confirm the theoretical results. The resulting adaptive mesh refinement procedures in 3d recover the adaptive convergence rates known for elliptic problems.

Introduction

The efficient numerical treatment of boundary integral equations using adaptive mesh refinement procedures has been extensively investigated for the numerical solution of homogeneous elliptic problems in unbounded domains [START_REF] Carstensen | Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasi-uniform meshes[END_REF][START_REF] Carstensen | A posteriori error estimates for boundary element methods[END_REF]. See [START_REF] Gwinner | Advanced Boundary Element Methods -Treatment of Boundary Value, Transmission and Contact Problems[END_REF] for a recent exposition.

In this article we investigate the extension of the a posteriori error analysis and adaptive mesh refinement procedures to initial-boundary value problems for the wave equation, formulated as boundary integral equations in the time-domain [START_REF] Costabel | Time-dependent problems with the boundary integral equation method[END_REF][START_REF] Sayas | Retarded Potentials and Time Domain Boundary Integral Equations: A Road Map[END_REF]. We prove a reliable a posteriori error estimate of residual type for a large class of conforming discretizations. It is efficient for a time-domain boundary element method on a globally quasi-uniform mesh. The error estimate defines an adaptive mesh refinement procedure, which recovers the convergence rates known for time-independent screen problems.

There has been recent interest in the solution of such problems on adapted meshes. Similar to the elliptic case, singularities of the solution may appear at singular points of the boundary, as discussed in [START_REF] Kokotov | The Neumann problem for the wave equation in a cone[END_REF][START_REF] Kokotov | Diffraction on a cone: The asymptotics of solutions near the vertex[END_REF][START_REF] Plamenevskiǐ | On the Dirichlet problem for the wave equation in a cylinder with edges[END_REF], and in trapping regions. Time-independent graded meshes have been shown to recover quasi-optimal convergence rates for edge and corner singularities [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF]. First steps towards time-adaptivity are due to Sauter and Veit [START_REF] Sauter | Adaptive Time Discretization for Retarded Potentials[END_REF] in 2 dimensions, and also convolution quadrature methods with graded, non-adaptively chosen time steps have been studied, for example in [START_REF] Sauter | Convolution quadrature for the wave equation with impedance boundary conditions[END_REF]. Gläfke [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF] showed first results towards space-time refinements in 2 dimensions, and in unpublished work Abboud uses ZZ error indicators for computations with space-adaptive mesh refinements for screen problems.

To describe the main results, we consider the wave equation

∂ 2
t u -∆u = 0 , u = 0 for t ≤ 0 , in the complement of a polyhedral domain or screen Ω ⊂ R d , with an emphasis on the challenging case d = 3. On the boundary Γ = ∂Ω both Dirichlet and acoustic boundary conditions, u = f , respectively ∂ ν u -α∂ t u = f , are considered. Here f is given, ν is the outer unit normal vector to Γ, and 0 < α, α -1 ∈ L ∞ (Γ). Following Bamberger and Ha Duong [START_REF] Bamberger | Formulation variationnelle espace-temps pour le calcul par potentiel retard de la diffraction d'une onde acoustique[END_REF], we recast the boundary problem as a time dependent boundary integral equation. The Dirichlet problem is equivalent to a hyperbolic variant of Symm's integral equation:

Vφ(t, x) = R + ×Γ
G(tτ, x, y)φ(τ, y) dτ ds y = f (t, x) .

(

) 1
Here φ is sought in a space-time anisotropic Sobolev space H 1 σ (R + , H -1 2 (Γ)), and G is a fundamental solution of the wave equation, G(ts, x, y) = H(ts -|x -y|)

2π (t -s) 2 + |x -y| 2 in 2d, (2)
G(ts, x, y) = δ(ts -|x -y|) 4π|x -y| in 3d.

(

) 3
Here H is the Heaviside function and δ the Dirac distribution. Our results apply, in particular, to a Galerkin discretization of the weak form of (1) in a subspace V ⊂ H 1 σ (R + , H -1 2 (Γ)),

Theorem A: Let φ ∈ H 1 σ (R + , H -1 2 (Γ)) be the solution to (4), and let φ h,∆t ∈ H

1 σ (R + , H -1 2 (Γ)) such that R = ∂ t f -V∂ t φ h,∆t ∈ H 0 σ (R + , H 1 (Γ)). Then φ -φ h,∆t 2 0,-1 2 ,Γ, * i,∆ max{∆t, h ∆ } R 2 0,1,[t i ,t i+1 )×∆ . (5)
Let Γ be closed and polyhedral. For a globally quasi-uniform mesh on Γ, let V = W 0 h,∆t the tensor product of cubic splines in time with piecewise constant functions in space. If φ h,∆t ∈ V is a Galerkin solution of (4) in V and φ ∈ H 2 σ (R + , H -1 2 (Γ)), then

max{∆t, h} R 2 0,1-ǫ,Γ φ -φ h,∆t 2 2,-1 2 ,Γ . (6)
The upper bound [START_REF] Carstensen | A posteriori error estimates for boundary element methods[END_REF] is obtained in Corollary 4.5, the lower bound [START_REF] Carstensen | A posteriori error estimate and h-adaptive algorithm on surfaces for Symm's integral equation[END_REF] in Theorems 5.1 and 6.1. Our numerical results illustrate the a posteriori error estimate of Theorem A for time-domain boundary elements based on [START_REF] Carstensen | Averaging techniques for the effective numerical solution of Symm's integral equation of the first kind[END_REF], but the estimate also applies, for example, to convolution quadrature methods [START_REF] Sayas | Retarded Potentials and Time Domain Boundary Integral Equations: A Road Map[END_REF].

Note the loss of time derivatives between the upper and lower bound of the error, in the first Sobolev index. The loss is well-known for error estimates for hyperbolic problems [START_REF] Duong | On retarded potential boundary integral equations and their discretizations[END_REF]. Our arguments generalize to give reliable a posteriori estimates for the acoustic boundary problem, see Section 3.2.

The residual error estimate from Theorem A is used to define adaptive mesh refinements in space, based on the four steps Solve, Estimate, Mark, Refine. Numerical experiments confirm the efficiency and reliability of the estimate in examples. For screen problems they recover the convergence rates known for elliptic problems. This work builds on the numerical analysis of adaptive boundary element methods for the Laplace equation, both for Symm's integral equation and the hypersingular equation [START_REF] Carstensen | Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasi-uniform meshes[END_REF][START_REF] Carstensen | A posteriori error estimates for boundary element methods[END_REF][START_REF] Carstensen | A posteriori error estimate and h-adaptive algorithm on surfaces for Symm's integral equation[END_REF][START_REF] Carstensen | Residual-based a posteriori error estimate for hypersingular equation on surfaces[END_REF]. Work on different types of error indicators in the time-independent case includes ZZ [START_REF] Carstensen | Averaging techniques for the effective numerical solution of Symm's integral equation of the first kind[END_REF] and Faermann indicators [START_REF] Faermann | Local a-posteriori error indicators for the Galerkin discretization of boundary integral equations[END_REF][START_REF] Faermann | Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods. II. The three-dimensional case[END_REF]. We also mention [START_REF] Wendland | Adaptive boundary element methods for strongly elliptic integral equations[END_REF] for an earlier approach. A comparison of different indicators in the time-domain will be the subject of future research.

Structure of this article: Section 2 recalls the boundary integral operators associated to the wave equation as well as their mapping properties between suitable space-time anisotropic Sobolev spaces. The Sobolev spaces are discretized using tensor products of piecewise polynomials in space and time. Section 3 presents a corresponding space-time discretization for the formulation of the Dirichlet problem in terms of the single layer operator and derives a reliable a posteriori error estimate for globally quasi-uniform meshes, using a canonical approach which will readily adapt to other settings. A second subsection analyzes an acoustic boundary problem, a system of equations involving in addition the double layer, adjoint double layer and hypersingular operators. Section 4 then localizes the space-time Sobolev norm to derive the general upper estimate for the Dirichlet problem in Theorem A. The upper estimates are complemented by a lower bound for the error of a Galerkin approximation on globally quasi-uniform meshes in Section 5. The final step of this proof is a lower bound for the best approximation in Section 6. Section 7 discusses some algorithmic properties of the implementation, before numerical experiments are used to confirm the theoretical results in Section 8.

Preliminaries and discretization

In addition to the single layer operator V, for acoustic boundary problems we require its normal derivative K ′ , the double layer operator K and hypersingular operator W for x ∈ Γ, t > 0:

Vϕ(t, x) = R + ×Γ G(t -τ, x, y) ϕ(τ, y) dτ ds y , (7)
Kϕ(t, x) = R + ×Γ ∂G ∂n y (t -τ, x, y) ϕ(τ, y) dτ ds y , (8)
K ′ ϕ(t, x) = R + ×Γ ∂G ∂n x (t -τ, x, y) ϕ(τ, y) dτ ds y , (9)
Wϕ(t, x) = R + ×Γ ∂ 2 G ∂n x ∂n y (t -τ, x, y) ϕ(τ, y) dτ ds y . (10)
Space-time anisotropic Sobolev spaces on the boundary Γ provide a convenient setting to study the mapping properties of layer potentials. See [START_REF] Duong | On retarded potential boundary integral equations and their discretizations[END_REF][START_REF] Gimperlein | A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space[END_REF] for a detailed exposition. To define them, if ∂Γ = ∅, first extend Γ to a closed, orientable Lipschitz manifold Γ.

On Γ one defines the usual Sobolev spaces of supported distributions:

H s (Γ) = {u ∈ H s ( Γ) : supp u ⊂ Γ} , s ∈ R . Furthermore, H s (Γ) is the quotient space H s ( Γ)/ H s ( Γ \ Γ).
To write down an explicit family of Sobolev norms, introduce a partition of unity α i subordinate to a covering of Γ by open sets B i . For diffeomorphisms ϕ i mapping each B i into the unit cube ⊂ R n , a family of Sobolev norms is induced from R n , with parameter ω ∈ C \ {0}:

||u|| s,ω, Γ = p i=1 R n (|ω| 2 + |ξ| 2 ) s |F (α i u) • ϕ -1 i (ξ)| 2 dξ 1 2
.

The norms for different ω ∈ C \ {0} are equivalent and F denotes the Fourier transform. They induce norms on H s (Γ), ||u|| s,ω,Γ = inf v∈ H s ( Γ\Γ) ||u + v|| s,ω, Γ and on H s (Γ), ||u|| s,ω,Γ, * = ||e + u|| s,ω, Γ . We write H s ω (Γ) for H s (Γ), respectively H s ω (Γ) for H s (Γ), when a norm with a specific ω is fixed. e + extends the distribution u by 0 from Γ to Γ. As the norm ||u|| s,ω,Γ, * corresponds to extension by zero, while ||u|| s,ω,Γ allows extension by an arbitrary v, ||u|| s,ω,Γ, * is stronger than ||u|| s,ω,Γ . Like in the time-independent case the norms are not equivalent whenever s ∈ 1 2 + Z [START_REF] Gwinner | Advanced Boundary Element Methods -Treatment of Boundary Value, Transmission and Contact Problems[END_REF]. We now define a class of space-time anisotropic Sobolev spaces: Definition 2.1. For σ > 0 and r, s ∈ R define .

H r σ (R + , H s (Γ)) = {u ∈ D ′ + (H s (Γ)) : e -σt u ∈ S ′ + (H s (Γ)) and ||u|| r,s,σ,Γ < ∞} , H r σ (R + , H s (Γ)) = {u ∈ D ′ + ( H s (Γ)) : e -σt u ∈ S ′ + ( H s (Γ))
They are Hilbert spaces, and we note that the basic case r = s = 0 is the weighted L 2 -space with scalar product u, v σ := ∞ 0 e -2σt Γ uvds x dt. Because Γ is Lipschitz, like in the case of standard Sobolev spaces these spaces are independent of the choice of α i and ϕ i when |s| ≤ 1. Using variational arguments, precise mapping properties are well-known for the layer potentials between Sobolev spaces related to the energy: Theorem 2.2 ( [START_REF] Gimperlein | A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space[END_REF]). The following operators are continuous for r ∈ R:

V : H r+1 σ (R + , H -1 2 (Γ)) → H r σ (R + , H 1 2 (Γ)) , K ′ : H r+1 σ (R + , H -1 2 (Γ)) → H r σ (R + , H -1 2 (Γ)) , K : H r+1 σ (R + , H 1 2 (Γ)) → H r σ (R + , H 1 2 (Γ)) , W : H r+1 σ (R + , H 1 2 (Γ)) → H r σ (R + , H -1 2 (Γ)) .
More generally, outside the energy spaces one has: [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF]). The following operators are continuous for r ∈ R, s ∈ (-1 2 , 1 2 ):

Theorem 2.3 ([
V : H r+1 σ (R + , H -1 2 +s (Γ)) → H r σ (R + , H 1 2 +s (Γ)) , K ′ : H r+2 σ (R + , H -1 2 +s (Γ)) → H r σ (R + , H -1 2 +s (Γ)) , K : H r+2 σ (R + , H 1 2 +s (Γ)) → H r σ (R + , H 1 2 +s (Γ)) , W : H r+3 σ (R + , H 1 2 +s (Γ)) → H r σ (R + , H -1 2 +s (Γ)) .
For simplicity, we will use here a hypersurface composed of N s triangular facets such that Γ = ∪ Ns i=1 T i , each element T i is closed with int(T i ) = ∅, and for distinct T i , T j ⊂ Γ the intersection int(T i ) ∩ int(T j ) = ∅.

For the time discretization we consider a uniform decomposition of the time interval R + into subintervals I n = [t n-1 , t n ) with time step |I n | = ∆t, such that t n = n∆t, n = 0, 1, . . . . We choose a basis {ϕ p j } of the space V p h of piecewise polynomial functions of degree p in space (continuous if p ≥ 1) and a basis {β j,q } of the space V q ∆t of piecewise polynomial functions of degree of q in time (continuous and vanishing at t = 0 if q ≥ 1). In addition to V q ∆t we also require the space W ∆t ⊂ V 3 ∆t of cubic splines. We denote the finite temporal mesh

{[0, t 1 ), [t 1 , t 2 ), • • • , [t Nt-1 , t Nt )} by T T , the spatial mesh {T 1 , • • • , T Ns } for Γ by T S .
We consider the tensor product of the approximation spaces in space and time, V p h and V q ∆t , associated to the space-time mesh T S,T = T S × T T , and we write

V p,q h,∆t = V p h ⊗ V q ∆t .
We consider the orthogonal projections Π ∆t from L 2 (R + ) to V q ∆t , resp. Π h from L 2 (Γ) to V p h . See [START_REF] Gimperlein | A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space[END_REF] for a discussion of their properties and those of their composition Π h,∆t . Furthermore, we define W p h,∆t = V p h ⊗ W ∆t . One obtains as in Proposition 3.54 of [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF]:

Lemma 2.4. Let f ∈ H r σ (R + , H m (Γ)∩ H s (Γ)), 0 < m ≤ q + 1, 0 < r ≤ p + 1, s ≤ r, |l| ≤ 1 2 such that ls ≥ 0. Then if l, s ≤ 0 f -Π h • Π ∆t f s,l,Γ ≤ C(h α + (∆t) β )||f || r,m,Γ , where α = min{m -l, m -m(l+s) m+r }, β = min{m + r -(l + s), m + r -m+r m l}. If l, s > 0, β = m + r -(l + s).
For the lower bounds, we further require an auxiliary projection P h,∆t from H

2 (R + , H -1 (Γ)) into the space of ansatz functions ⊆ H 2 (R + , L 2 (Γ)) such that P h,∆t is H 2 (R + , L 2 (Γ))-stable.
For the construction of P h,∆t = P ∆t Q * h we use a projection P ∆t on the space of cubic splines W ∆t ⊂ V 3 ∆t in time and a Galerkin projection Q * h in space. More precisely, we consider the projection

Q h : L 2 (Γ) → V p h defined as follows: For φ ∈ L 2 (Γ), we let Q h φ ∈ V p
h be the unique solution to the variational problem

Q h φ, ψ h = φ, ψ h for all ψ h ∈ V p h . Q h is bounded on L 2 (Γ) and self-adjoint, Q h = Q * h . Therefore, if Q h restricts to a bounded operator on H s ω (Γ), Q h = Q * h extends
Q h φ s,ω,Γ ≤ C φ s,ω,Γ .
For all φ ∈ H -s ω (Γ), we have:

Q h φ -s,ω,Γ, * ≤ C φ -s,ω,Γ, * .
Proof. For the first assertion we note that it is clear for s = 0. We show it for s = 1, and the case for general s ∈ [0, 1] follows from interpolation.

Let Q h,1 : H 1 ω (Γ) → V p h ⊂ H 1 ω (Γ) defined by Q h,1 φ, ψ h 1,ω,Γ = φ, ψ h 1,ω,Γ for all ψ h ∈ V p h . By definition, Q h,1 is bounded on H 1 ω (Γ). We use the Aubin-Nitsche trick to show that φ -Q h,1 φ L 2 (Γ) h φ H 1 ω (Γ) . More precisely, if g ∈ L 2 (Γ) denote by φ g ∈ H 1 ω (Γ) the solution to the elliptic equation ψ, φ g 1,ω,Γ = g, ψ L 2 (Γ) , for all ψ ∈ H 1 ω (Γ) . By elliptic regularity φ g ∈ H 2 ω (Γ), with φ g 2,ω,Γ g L 2 (Γ) , as Γ does not have reentrant corners. We observe φ -Q h,1 φ L 2 (Γ) = sup 0 =g∈L 2 (Γ) g, φ -Q h,1 φ L 2 (Γ) g L 2 (Γ) = sup 0 =g∈L 2 (Γ) φ -Q h,1 φ, φ g 1,ω,Γ g L 2 (Γ) = sup 0 =g∈L 2 (Γ) inf 0 =ψ h ∈V p h φ -Q h,1 φ, φ g -ψ h 1,ω,Γ g L 2 (Γ) ≤ φ -Q h,1 φ 1,ω,Γ sup 0 =g∈L 2 (Γ) inf 0 =ψ h ∈V p h φ g -ψ h 1,ω,Γ g L 2 (Γ) . Because φ g ∈ H 2 ω (Γ), inf 0 =ψ h ∈V p h φ g -ψ h 1,ω,Γ h ∇ Γ φ g 1,ω,Γ h g L 2 (Γ) and φ -Q h,1 φ 1,ω,Γ φ 1,ω,Γ , we conclude φ -Q h,1 φ L 2 (Γ) h φ 1,ω,Γ .
Combining this with the triangle inequality, the inverse inequality and the already shown assertion when s = 0, we obtain the first assertion:

Q h φ 1,ω,Γ ≤ Q h,1 φ 1,ω,Γ + Q h φ -Q h,1 φ 1,ω,Γ φ 1,ω,Γ + h -1 Q h φ -Q h,1 φ L 2 (Γ) = φ 1,ω,Γ + h -1 Q h (φ -Q h,1 φ) L 2 (Γ) φ 1,ω,Γ + h -1 φ -Q h,1 φ L 2 (Γ) φ 1,ω,Γ .
The second assertion follows from the first one, with the same constant C. Indeed,

Q h φ -s,ω,Γ, * = sup 0 =ψ∈H s ω (Γ) Q h φ, ψ L 2 (Γ) ψ s,ω,Γ = sup 0 =ψ∈H s ω (Γ) φ, Q h ψ L 2 (Γ) ψ s,ω,Γ ≤ C φ -s,ω,Γ, * .
Using the Laplace transform to transfer the result into the time-domain, we obtain

Q h φ r,-s,Γ, * ≤ C φ r,-s,Γ, *
for s ∈ [0, 1] and all r. Further note that the interpolation operator P ∆t in time for cubic splines is continuous on H 2 functions. The continuity extends to H 2 (R + , H -s (Γ)). We conclude:

Lemma 2.6. Consider a quasi-uniform mesh, s ∈ [0, 1] and P h,∆t = P ∆t Q h . For all φ ∈ H 2 (R + , H -s (Γ)), we have

P h,∆t φ 2,-s,Γ, * ≤ C φ 2,-s,Γ, * .
The approximation properties of P h,∆t analogous to Lemma 2.4 follow as in [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF], Proposition 3.54.

A posteriori error estimates -reliability 3.1 Dirichlet problem

We recall the basic properties of the bilinear form

B D (φ, ψ) = R + Γ V∂ t φ(t, x) ψ(t, x) ds x d σ t of the Dirichlet problem.
As shown in [START_REF] Duong | On retarded potential boundary integral equations and their discretizations[END_REF], the bilinear form is continuous, and also weakly coercive:

Proposition 3.1. For every φ, ψ ∈ H 1 σ (R + , H -1 2 (Γ)
) there holds:

|B D (φ, ψ)| φ 1,-1 2 ,Γ, * ψ 1,-1 2 ,Γ, * and φ 2 0,-1 2 ,Γ, * B D (φ, φ).
We consider a conforming Galerkin discretization of the Dirichlet problem (4) in a subspace V ⊂ H 1 σ (R + , H -1 2 (Γ)), which reads as follows: Find φ h,∆t ∈ V such that

B D (φ h,∆t , ψ h,∆t ) = ∂ t f, ψ h,∆t , (11)
for all ψ h,∆t ∈ V .

The well-posedness of the continuous and discretized problems are a basic consequence of Proposition 3.1:

Corollary 3.2. Let f ∈ H 2 σ (R + , H 1 2 (Γ)).
Then the Dirichlet problem (4) and its discretization [START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF]

admit unique solutions φ ∈ H 1 σ (R + , H -1 2 (Γ)), φ h,∆t ∈ V , and the estimates φ 1,-1 2 ,Γ, * , φ h,∆t 1,-1 2 ,Γ, * f 2, 1 2 ,Γ hold.
We note the Galerkin orthogonality:

B D (φ -φ h,∆t , ψ h,∆t ) = 0 ∀ψ h,∆t ∈ V .
Using ideas going back to Carstensen [START_REF] Carstensen | Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasi-uniform meshes[END_REF] and Carstensen and Stephan [START_REF] Carstensen | A posteriori error estimates for boundary element methods[END_REF] for the boundary element method for elliptic problems, we obtain an a posteriori error estimate for the Galerkin solution to the Dirichlet problem on globally quasi-uniform meshes.

Theorem 3.3. Given V = V p,q h,∆t associated to a globally quasi-uniform discretization of Γ, let φ ∈ H 1 σ (R + , H -1 2 (Γ)), φ h,∆t ∈ V the solutions to (4), resp. (11). Assume that R = ∂ t f -V∂ t φ h,∆t ∈ H 0 σ (R + , H 1 (Γ)). Then φ -φ h,∆t 2 0,-1 2 ,Γ, * R 0,1,Γ ∆t ∂ t R 0,0,Γ + h • ∇R 0,0,Γ max{∆t, h}( ∂ t R 2 0,0,Γ + ∇R 2 0,0,Γ ) .
Proof of Theorem 3.3. We first note that for all ψ h,∆t ∈ V p,q h,∆t

φ -φ h,∆t 2 0,-1 2 ,Γ, * B D (φ -φ h,∆t , φ -φ h,∆t ) = R + Γ ∂ t f (φ -φ h,∆t ) ds x d σ t -B D (φ h,∆t , φ -φ h,∆t ) = R + Γ ∂ t f (φ -ψ h,∆t ) ds x d σ t -B D (φ h,∆t , φ -ψ h,∆t ) = R + Γ (∂ t f -V∂ t φ h,∆t )(φ -ψ h,∆t ) ds x d σ t .
The last term may be estimated by:

R + Γ (∂ t f -V φh,∆t )(φ -ψ h,∆t ) ds x d σ t ≤ R 0, 1 2 ,Γ φ -ψ h,∆t 0,- 1 2 ,
Γ, * . We use ψ h,∆t = φ h,∆t together with the interpolation inequality

R 2 0, 1 2 ,Γ ≤ R 0,0,Γ R 0,1,Γ .
As the residual is perpendicular to V p,q h,∆t ,

R 2 0,0,Γ = R, R = R, R -ψ h,∆t ≤ R 0,0,Γ R -ψ h,∆t 0,0,Γ
for all ψ h,∆t ∈ V p,q h,∆t , we obtain

R 0,0,Γ ≤ inf{ R -ψ h,∆t 0,0,Γ : ψ h,∆t ∈ V p,q h,∆t } .
Choosing ψ h,∆t = Π h,∆t R, based on the interpolation operator defined earlier, we obtain

R 0,0,Γ ∆t ∂ t R 0,0,Γ + h • ∇R 0,0,Γ .
The theorem follows.

Acoustic boundary problems

Recall the wave equation with inhomogeneous acoustic boundary conditions

∂ 2 t u -∆u = 0 on R + × Ω , ∂ ν u -α∂ t u = f on R + × Γ , u = 0 for t ≤ 0 .
For scattering problems f = -∂ ν u inc + α∂ t u inc is determined from an incoming wave u inc . For a finite or infinite time interval [0, T ] we introduce the bilinear form

a T ((ϕ, p), (ψ, q)) = T 0 Γ α φ ψ + 1 α pq + K ′ p ψ -Wϕ ψ + V ṗq -K φq ds x dt . (12)
With

l(ψ, q) = T 0 Γ F ψ ds x dt + ∞ 0 Γ Gq α ds x dt , (13)
where F = -2∂ ν u inc , G = 2α∂ t u inc , we consider the variational formulation for the wave equation in R 3 with acoustic boundary conditions on Γ:

Find (ϕ, p) ∈ H 1 ([0, T ], H 1 2 (Γ)) × H 1 ([0, T ], L 2 (Γ)) such that a T ((ϕ, p), (ψ, q)) = l(ψ, q) (14) for all (ψ, q) ∈ H 1 ([0, T ], H 1 2 (Γ)) × H 1 ([0, T ], L 2 (Γ)).
Note that σ may be set to 0 in the definition of the Sobolev spaces when T < ∞. The acoustic problem is equivalent to the wave equation with acoustic boundary conditions [START_REF] Duong | On retarded potential boundary integral equations and their discretizations[END_REF].

Its discretization reads: Find

(ϕ h,∆t , p h,∆t ) ∈ V p 1 ,q 1 h,∆t × V p 2 ,q 2 h,∆t such that a T ((ϕ h,∆t , p h,∆t ), (ψ h,∆t , q h,∆t )) = l(ψ h,∆t , q h,∆t ) (15)
for all (ψ h,∆t , q h,∆t ) ∈ V p 1 ,q 1 h,∆t × V p 2 ,q 2 h,∆t .

The following well-posedness holds:

Proposition 3.4. Let F ∈ H 2 ([0, T ], H -1 2 (Γ)), G ∈ H 1 ([0, T ], H 0 (Γ)).
Then the weak form [START_REF] Gimperlein | A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space[END_REF] of the acoustic problem and its discretization [START_REF] Gimperlein | Time domain boundary element methods for the Neumann problem: Error estimates and acoustic problems[END_REF]

admit unique solutions (ϕ, p) ∈ H 1 ([0, T ], H 1 2 (Γ)) × H 1 ([0, T ], L 2 (Γ)), resp. (ϕ h,∆t , p h,∆t ) ∈ V p 1 ,q 1 h,∆t × V p 2 ,q 2
h,∆t , which depend continuously on the data.

We specifically note that the bilinear form a T satisfies a (weaker) coercivity estimate:

p 2 0,0,Γ + φ 2 0,0,Γ a T ((ϕ, p), (ϕ, p)) .
This follows from Equation (64) of [START_REF] Duong | On retarded potential boundary integral equations and their discretizations[END_REF],

a T ((ϕ, p), (ϕ, p)) = 2E(T ) + T 0 Γ α(∂ t ϕ)(∂ t ϕ) + 1 α pq ds x dt ,
where 2 dx is the total energy at time T . We state a simple a posteriori estimate.

E(T ) = 1 2 R d \Ω (∂ t u) 2 + (∇u)
Theorem 3.5. Let (ϕ, p) ∈ H 1 ([0, T ], H 1 2 (Γ)) × H 1 ([0, T ], L 2 (Γ)
) be the solution to [START_REF] Gimperlein | A priori error estimates for a time-dependent boundary element method for the acoustic wave equation in a half-space[END_REF], (ϕ h,∆t , p h,∆t ) ∈ V p 1 ,q 1 h,∆t × V p 2 ,q 2 h,∆t the solution to the discretization [START_REF] Gimperlein | Time domain boundary element methods for the Neumann problem: Error estimates and acoustic problems[END_REF]. Assume that

R 1 = F -α φh,∆t + 2K ′ p h,∆t -2Wϕ h,∆t ∈ L 2 ([0, T ], L 2 (Γ)) , R 2 = G + α -1 p h,∆t + 2V ṗh,∆t -2K φh,∆t ∈ L 2 ([0, T ], L 2 (Γ)) .
Then pp h,∆t 0,0,Γ + φφh,∆t 0,0,Γ R 1 0,0,Γ + R 2 0,0,Γ .

Proof. For every (ψ h,∆t , q h,∆t ) ∈ V p 1 ,q 1 h,∆t × V p 2 ,q 2 h,∆t we have

p -p h,∆t 2
0,0,Γ + φ -φh,∆t 2 0,0,Γ a T ((ϕ -ϕ h,∆t , p -p h,∆t ), (ϕ -ϕ h,∆t , p -p h,∆t )) = R 1 , φ -ψh,∆t + R 2 , p -q h,∆t ≤ R 1 0,0,Γ φ -ψh,∆t 0,0,Γ + R 2 0,0,Γ p -q h,∆t 0,0,Γ ≤ ( R 1 0,0,Γ + R 2 0,0,Γ )× ( p -q h,∆t 0,0,Γ + φ -ψh,∆t 0,0,Γ ) .
The assertion is obtained by choosing (ψ h,∆t , q h,∆t ) = (ϕ h,∆t , p h,∆t ).

Naturally, for a quasi-uniform discretization of Γ and under stronger assumptions on R 1 , R 2 we may obtain powers of h and ∆t on the right hand side by the following argument:

As in the proof of Theorem 3.3, R 1 , ˙ ψ h,∆t = R 2 , q h,∆t = 0 for all ( ψ h,∆t , q h,∆t ) ∈ V p 1 ,q 1 h,∆t × V p 2 ,q 2 h,∆t . Hence R 2 2 0,0,Γ = R 2 , R 2 = R 2 , R 2 -q h,∆t ≤ R 2 0,0,Γ R 2 -q h,∆t 0,0,Γ . Choosing q h,∆t = Π h,∆t R 2 yields R 2 0,0,Γ ∆t ∂ t R 2 0,0,Γ + h • ∇ Γ R 2 0,0,Γ + ∆t h • ∇ Γ ∂ t R 2 0,0,Γ provided R 2 ∈ H 1 ([0, T ], H 1 (Γ)). Assuming R 1 ∈ H 1 ([0, T ], H 1 (Γ)), we similarly have R 1 2 0,0,Γ = R 1 , R 1 = R 1 , R 1 - ˙ ψ h,∆t ≤ R 1 1 2 -s,0,Γ t 0 R 1 -ψ h,∆t 1 2 +s,0,Γ . Choosing ψ h,∆t = Π h,∆t t 0 R 1 and s = 1 2 results in R 1 0,0,Γ ∆t ∂ t R 1 0,0,Γ + h • ∇ Γ R 1 0,0,Γ + ∆t h • ∇ Γ ∂ t R 1 0,0,Γ .
Altogether, pp h,∆t 0,0,Γ + φφh,∆t 0,0,Γ + ϕϕ h,∆t 0,

1 2 ,Γ R 1 0,0,Γ + R 2 0,0,Γ 2 i=1 ∆t ∂ t R i 0,0,Γ + h • ∇R i 0,0,Γ + ∆t h • ∇∂ t R i 0,0,Γ .

Error estimates for general discretizations

This section generalizes the results for the single layer potential without any assumptions on the underlying meshes. We recall Lemma 3 in [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF]:

Lemma 4.1. Let f 1 , . . . , f n ∈ H r σ (R + , H s (Γ)), 0 ≤ s ≤ 1, r ≥ 0, such that f j f k = 0 for j = k.
Let ω j be the interior of the support of f j with ω j = supp f j . Then

n j=1 f j 2 r,s,Γ, * ≤ n j=1 f j 2 r,s,ω j , * .
The lemma will be applied to a finite partition of unity Φ given by non-negative Lipschitz functions {φ j } M j=1 on Γ such that M j=1 φ j = 1.

Definition 4.2. The overlap of the partition of unity Φ is defined as K(Φ) = max j card{k :

φ k φ j = 0}.
For a partition of unity associated to a triangulation, M tends to infinity as the mesh size decreases, while the overlap may be much smaller. We note a crucial observation from [START_REF] Carstensen | A posteriori error estimate and h-adaptive algorithm on surfaces for Symm's integral equation[END_REF], Lemma 3.1: Lemma 4.3. Let Φ be a finite partition of unity of ∂Ω with overlap K(Φ). Then there exists a partition of {1, . . . , M } into K ≤ K(Φ) non-empty subsets M 1 , . . . , M K , such that

K j=1 M j = {1, . . . , M }, M j ∩ M k = ∅ if j = k
and for all l ∈ {1, . . . , K} and j, k ∈ M l with j = k, φ j φ k = 0 on [0, T ] × ∂Ω.

Similar to the partition of unity in space, associated to the decomposition of the time axis, R + = ∞ l=0 [l∆t, (l+1)∆t), we consider a partition of unity {ψ l } ∞ l=0 on R + , ∞ l=0 ψ l = 1. We may choose it so that ψ l is supported in the interval ((l -1

2 )∆t, (l + 3 2 )∆t), and therefore

ψ l ψ l ′ = 0 whenever l ∈ M 0 = 2N, l ′ ∈ M 1 = 2N + 1.
We obtain a partition of unity { ψ l,j = ψ l ⊗ φ j } l,j on R + × ∂Ω.

Theorem 4.4. Let Γ ′ ⊂ Γ be connected. and let Φ be a finite partition of unity with overlap K(Φ). Then for any f ∈ H r σ (R + , H s (Γ ′ )) and any 0 ≤ s ≤ 1, we have

f 2 r,s,Γ ′ , * ≤ 2K(Φ) ∞ l=0 M j=1 f ψ l,j 2(1-s) r,0,Γ ′ f ψ l,j 2s r,1,Γ ′ .
Proof. We show

f 2 r,s,Γ ′ , * ≤ 2K(Φ) ∞ l=0 M j=1 f ψ l,j 2 r,s,Γ ′ , * . (16)
The assertion then follows from the interpolation estimate f ψ l,j r,s,Γ ′ , * f ψ l,j 1-s r,0,Γ ′ f ψ l,j s r,1,Γ ′ .

To show [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF], we consider a partition M 1 , . . . , M K as in Lemma 4.3. Then f = ∞ l=0 M j=1 ψ l,j f , so that

f 2 r,s,Γ ′ , * = 1 m=0 l∈ Mm K k=1 j∈M k ψ l,j f 2 r,s,Γ ′ , * ≤ 2K 1 m=0 K k=1 l∈ Mm j∈M k ψ l,j f 2 r,s,Γ ′ , * .
With Lemma 4.1 and I l = ((l -1)∆t, (l + 2)∆t),

l∈ Mm j∈M k ψ l,j f 2 r,s,Γ ′ , * ≤ j∈M k l∈ Mm ψ l,j f 2 r,s,ω j , * ≤ l∈ Mm j∈M k ψ l,j f 2 r,s,I l ×ω j , * , so that f 2 r,s,Γ ′ , * ≤ 2K 1 m=0 l∈ Mm K k=1 j∈M k ψ l,j f 2 r,s,I l ×ω j , * = 2K ∞ l=0 M j=1
ψ l,j f 2 r,s,I l ×ω j , * , which shows [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF].

Together with Friedrichs inequality

f ψ l,j r,1,Γ ′ σ ∇(f ψ l,j ) r,0,Γ ′ + ∂ t (f ψ l,j ) r,0,Γ ′ one obtains f ψ l,j 1-s r,0,Γ ′ f ψ l,j s r,1,Γ ′ σ d 2(1-s) j (1 + d 2 j ) s ∇(f ψ l,j ) 2 r,0,Γ ′ + ∂ t (f ψ l,j ) 2 r,0,Γ ′ .
Here d j = max{δ j , ∆t}, where δ j is the width of the support of φ j , defined as the smallest number such that the following is true: There exists a direction n ∈ R 3 , |n| = 1, such that for all x ∈ R 3 and for each plane H perpendicular to n with x ∈ H, the intersection supp φ j ∩ H is a Lipschitz curve of length ≤ δ j . We now use V∂ t (φφ h,∆t ) = R, the coercivity of V∂ t in Proposition 3.1, and Theorem 4.4 with s = 1 2 . The following a posteriori error estimate follows:

Corollary 4.5. Let φ ∈ H 1 σ (R + , H -1 2 (Γ)
) be the solution to (4), and let

φ h,∆t ∈ H 1 σ (R + , H -1 2 (Γ)) such that R = ∂ t f -V∂ t φ h,∆t ∈ H 0 σ (R + , H 1 (Γ)). Then φ-φ h,∆t 2 0,-1 2 ,Γ, * σ j δ j ∇R 2 0,0,I l ×∆ j + ∂ t R 2 0,0,I l ×∆ j ≃ l,∆ j max{∆t, h ∆ } R 2 0,1,I l ×∆ j .

Lower bounds

As for time-independent problems, for the discussion of lower bounds we restrict ourselves to globally quasi-uniform meshes.

Because of the different norms in the upper and lower bounds for B D in Proposition 3.1, the a posteriori estimate only satisfies a weak variant of efficiency: Provided φ ∈ H 2 σ (R + , H 0 (Γ)) we have for ε ∈ (0, 1):

max{∆t, h} -1-ε 2 φ -φ h,∆t 0,-1 2 ,Γ, * R 0,1-ε,Γ = V( φ -φh,∆t ) 0,1-ε,Γ φ -φ h,∆t 2,-ε,Γ ≤ φ -φ h,∆t 2,0,Γ .
A proof of the sharp estimate, ε = 0, would require sharp mapping properties of the layer potentials outside the energy spaces and will be pursued elsewhere.

As in the elliptic case, we aim to use the mapping properties of V together with approximation properties of the finite element spaces to recover the same spatial Sobolev index - 1 2 in the upper and lower estimates.

Theorem 5.1. Assume that R ∈ H 0 σ (R + , H 1 (Γ)) and that the ansatz functions V = W p h,∆t ⊂ H 2 σ (R + , H 0 (Γ)) ∩ V p,q h,∆t satisfy inf ψ h∆t ∈V φ -ψ h∆t 2,0,Γ ≃ max{h, ∆t} β (17)
for some β > 0. Then for all ε ∈ (0, 1)

R 0,1-ε,Γ max{h -1 2 , (∆t) -1 2 } φ -φ h∆t 2,-1 2 ,Γ .
Remark 5.2. In Theorem 6.1 below, the hypothesis (17) is verified using the singular expansion of the solution φ at the edges and corners.

Proof of Theorem 5.1. The theorem partly follows the functional analytic approach of [START_REF] Carstensen | Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasi-uniform meshes[END_REF] for the time-independent case. We consider the following quantities: the error of the best approximation Π h,∆t φ of φ in

H 2 σ (R + , L 2 (Γ)), E(φ, h, ∆t) = φ -Π h,∆t φ 2,0,Γ ,
where E(φ, h, ∆t) ≃ max{h, ∆t} β by [START_REF] Gwinner | Advanced Boundary Element Methods -Treatment of Boundary Value, Transmission and Contact Problems[END_REF]; the relative error of best approximation,

F (φ, h, ∆t) = φ -Π h,∆t φ 2,-1,Γ φ -Π h,∆t φ 2,0,Γ ;
and the inverse inequality [START_REF] Glaefke | Adaptive Methods for Time Domain Boundary Integral Equations[END_REF] G(h, ∆t) = sup

ψ h∆t =0 ψ h∆t 2,0,Γ ψ h∆t 2,-1,Γ max{h -1 , ∆t -1 }.
Their analysis relies on the projection operator P h,∆t from Lemma 2.6. Note that F (φ, h, ∆t) max{h, ∆t}. Indeed because φ ∈ H 2 σ (R + , H 0 (Γ)), by duality and the approximation properties of Π h,∆t :

φ -Π h,∆t φ 2,-1,Γ = sup η φ -Π h,∆t φ, η -Π h,∆t η η -2,1,Γ max{h, ∆t} φ -Π h,∆t φ 2,0,Γ .
This uses the estimate η -Π h,∆t η -2,0,Γ max{h, ∆t} η -2,1,Γ .

From the beginning of this section, we recall

R 0,1-ε,Γ φ -φ h,∆t 2,0,Γ ≤ φ -Π H,∆T φ 2,0,Γ + Π H,∆T φ -φ h,∆t 2,0,Γ .
Note that the first term is φ -Π H,∆T φ 2,0,Γ = E(φ, H, ∆T ). For the second term we observe that

Π H,∆T φ -φ h,∆t 2,0,Γ = P H,∆T (Π H,∆T φ -φ h,∆t ) 2,0,Γ G(H, ∆T ) Π H,∆T φ -φ h,∆t 2,-1,Γ , so that Π H,∆T φ -φ h,∆t 2,0,Γ G(H, ∆T ) 1/2 Π H,∆T φ -φ h,∆t 2,-1 2 ,Γ . We conclude R 0,1-ε,Γ E(φ, H, ∆T ) + G(H, ∆T ) 1/2 Π H,∆T φ -φ h,∆t 2,-1 2 ,Γ . Further, Π H,∆T φ -φ h,∆t 2,-1 2 ,Γ ≤ Π H,∆T φ -φ 2,-1 2 ,Γ + φ -φ h,∆t 2,-1 2 ,Γ . From the definition of F and interpolation, Π H,∆T φ -φ 2,-1 2 ,Γ is bounded by F (φ, H, ∆T ) 1/2 Π H,∆T φ -φ 2,0,Γ = F (φ, H, ∆T ) 1/2 E(φ, H, ∆T ) . To sum up, φ -φ h,∆t 2,0,Γ E(φ, H, ∆T ) + G(H, ∆T ) 1/2 F (φ, H, ∆T ) 1/2 E(φ, h, ∆t) + G(H, ∆T ) 1/2 φ -φ h,∆t 2,-1 2 ,Γ ≤ E(φ, H, ∆T ) E(φ, h, ∆t) φ -φ h,∆t 2,0,Γ (1 + G(H, ∆T ) 1/2 F (φ, H, ∆T ) 1/2 ) + G(H, ∆T ) 1/2 φ -φ h,∆t 2,-1 2 ,Γ , or, with δ = E(φ,H,∆T ) E(φ,h,∆t) (1 + F (φ, H, ∆T ) 1/2 G(φ, H, ∆T ) 1/2 ) and a constant C, φ -φ h,∆t 2,0,Γ G(H, ∆T ) 1/2 1 -Cδ φ -φ h,∆t 2,-1 2 ,Γ . If δ = E(φ,H,∆T ) E(φ,h,∆t) (1 + F (φ, H, ∆T ) 1/2 G(φ, H, ∆T ) 1/2 ) < 1 2C , we obtain R 0,1-ε,Γ G(H, ∆T ) 1/2 φ -φ h∆t 2,-1 2 ,Γ .
Here H and ∆T are sufficiently small compared to h and ∆t, and it remains to choose them so that δ < 1 2C . Set H = ρh and ∆T = ρ∆t for ρ ∈ (0, 1). Using that

F (φ, ρh, ρ∆t)G(φ, ρh, ρ∆t) ≃ max{h, ∆t} min{h, ∆t}
is uniformly bounded in ρ ∈ (0, 1), so that it suffices to show that E(φ,ρh,ρ∆t) E(φ,h,∆t) → 0 as ρ tends to 0. This follows from [START_REF] Gwinner | Advanced Boundary Element Methods -Treatment of Boundary Value, Transmission and Contact Problems[END_REF].

Best approximation and lower bounds

In this section verify the hypothesis (17) in Theorem 5.1 for polyhedral domains, by proving upper and lower bounds for the best approximation of the solution to the wave equation with Dirichlet boundary conditions.

Let Ω be a polyhedral domain and u a solution to the wave equation in Ω:

∂ 2 t u(t, x) -∆u(t, x) = 0 in R + t × Ω x , (18)
u(t, x) = g(t, x) on Γ = ∂Ω , (19)
u(0, x) = ∂ t u(0, x) = 0 in Ω. ( 20
)
The function u exhibits well-known singularities at non-smooth boundary points of the domain. Locally near an edge or a corner, Ω is of the form R + × K, where the base K ⊂ S 2 is a smooth or polygonal subset of the sphere. The solution may be decomposed into a leading part given by explicit singular functions plus less singular terms [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF][START_REF] Kokotov | The Neumann problem for the wave equation in a cone[END_REF][START_REF] Kokotov | Diffraction on a cone: The asymptotics of solutions near the vertex[END_REF][START_REF] Plamenevskiǐ | On the Dirichlet problem for the wave equation in a cylinder with edges[END_REF]. We refer to [START_REF] Kokotov | The Neumann problem for the wave equation in a cone[END_REF]Theorem 7.4 and Remark 7.5] for details in the case of the Neumann problem in a wedge, respectively [START_REF] Kokotov | Diffraction on a cone: The asymptotics of solutions near the vertex[END_REF]Theorem 4.1] for the Dirichlet problem in a cone and state the decomposition in terms of polar coordinates (r, θ) centered at the vertex (0, 0, 0):

u(t, x) = u 0 (t, r, θ) + χ(r)r λ a(t, θ) + χ(θ)b 1 (t, r)(sin(θ)) ν + χ( π 2 -θ)b 2 (t, r)(cos(θ)) ν , (21)
∂ n u(t, x) = ψ 0 (t, r, θ) + χ(r)r λ-1 a(t, θ) + χ(θ)b 1 (t, r)r -1 (sin(θ)) -ν + χ( π 2 -θ)b 2 (t, r)r -1 (cos(θ)) -ν . ( 22
)
Here, ν = π α , where α is the opening angle of the wedge, and λ = -1 2 + 1 4 + µ, where µ is the smallest eigenvalue of the Laplace-Beltrami operator with Dirichlet boundary conditions in the subdomain K of the sphere. χ, χ are cut-off functions and a, b j sufficiently regular. For generic problems the functions a, b 1 , b 2 are not identically zero.

From the representation formula, the wave equation translates into the boundary integral equation V φ = f , with f = (1 -K)g and solution φ = ∂ n u| Γ .

The main theorem concerning the approximation of φ is: Theorem 6.1. Assume that the coefficient functions a, b 1 , b 2 are not identically 0. Then E(φ, h, ∆t) ≃ max{h, ∆t} max{ν-1 2 ,λ} .

In particular, hypothesis (17) is satisfied. A similar result in the elliptic case was known if Γ is a curve, i.e. in dimension 2 [START_REF] Carstensen | Efficiency of a posteriori BEM-error estimates for first-kind integral equations on quasi-uniform meshes[END_REF].

The key step in the proof of Theorem 6.1 is to show the result for bilinear basis functions on a rectangular mesh. For simplicity of notation, we restrict to one boundary face of Γ and assume it is given by Q × {0} with corner of the domain at (0, 0, 0), where

Q = (0, 1) 2 = kl R kl , R kl = [x k-1 , x k ) × [x l-1 , x l ) and x k = kh.
Recall the following estimate from [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF]:

Lemma 6.2. Let -1 ≤ s ≤ 0, 0 ≤ r ≤ ρ ≤ p + 1, R = [0, h 1 ]×[0, h 2 ], u ∈ H ρ σ ([0, ∆t], H 1 (R)), Π p
t u the orthogonal projection onto piecewise polynomials in t of order p, and Π 0 x,y the orthogonal projection onto piecewise constant polynomials in space, Π 0

x,y u = 1

h 1 h 2 R
u(t, x, y)dy dx.

Then for U = Π p t Π 0 x,y u we have

u -U r,s,R, * (∆t) ρ-r max{h 1 , h 2 , ∆t} -s ∂ ρ t u L 2 ([0,∆t]×R) (23)
+ max{h 1 , h 2 , ∆t} -s h 1 u x L 2 ([0,∆t]×R) + h 2 u y L 2 ([0,∆t]×R) . If u(t, x, y) = u 1 (t, x)u 2 (y), u 1 ∈ H ρ σ ([0, ∆t], H 1 ([0, h 1 ])), u 2 ∈ H 1 ([0, h 2 ]) then u -U r,s,R, * (∆t) ρ-r max{h 1 , ∆t} -s ∂ ρ t u L 2 ([0,∆t]×R) + h 1-s 1 u x L 2 ([0,∆t]×R) + h 1-s 2 u y L 2 ([0,∆t]×R) .
Proof of 6.1. We use the upper bound E(φ, h, ∆t) max{h, ∆t} max{ν-1 2 ,λ} and first approximate the corner singularity f = r λ-1 a(t, θ). Here a is a function in

H ρ σ (R + , H 1 ([0, π/2])): With p| R kl = Π p t Π 0
x,y f and Lemma 6.2 one obtains

f -p 2 0,0,Q N k,l=1 k+l =2 (∆t) 2ρ ∂ ρ t f 2 0,0,R + h 2 k f x 2 0,0,R kl + h 2 l f y 2 0,0,R kl + f -p 2 0,0,R 11 . (24)
For k 2, l 2 the following estimate holds, with the zeta function ζ(s):

N k,l=2 h 2 k f x 2 0,0,R kl = N k=2 h 2 k ∞ 0 π/2 φ=0 w(t, θ) 2 dθ d σ t √ 2(k+1)h √ 2kh r 2λ-4 rdr ∞ k=2 h 2λ k 2λ-3 = ch 2λ ζ(3 -2λ) .
Here we have used

|f x (t, x, y)| r λ-2 w(t, θ) with w ∈ H ρ σ (R + , H 0 ([0, π/2])).
For k = 1, l > 1 (and analogously for k > 1, l = 1), in [START_REF] Petersdorff | Decompositions in edge and corner singularities for the solution of the Dirichlet problem of the Laplacian in a polyhedron[END_REF] we compute that

N l=2 h 2 1 f x 2 0,0,R 1l h 2 1 ∞ 0 h 1 x=0 1 y=0 |f x (t, x, y)| 2 dy dx d σ t h 2 1 ∞ 0 √ 2 r=h π/2 φ=0 a 2 (t, θ)dθ r 2λ-3 dr d σ t h 2 1 + h 2λ 1 .
Finally, in the corner k = 1, l = 1 the singular function

f ∈ H ρ σ (R + , H 0 (R 11 )), because λ > 0. Now f -p 0,0,R 11 f 0,0,R 11 ≃ h λ 1 .
We now consider the approximation of the edge singularities, which are of two types, fixing Q = [0, 1] 2 with singularity at the x-axis:

(i) f (t, x, y) = χ(x)b(t, x)y ν-1 with edge intensity factor b ∈ H ρ σ (R + , H 1 0 (R + )), (ii) f (t, x, y) = χ(x)b(t)x λ-ν y ν-1
with a corner singularity in the edge intensity factor.

We have

f -p 2 0,0,Q N k,l=1 f -p 2 0,0,R kl . (25)
First consider case (ii), where the time dependence factors out. Define 2 .

p 2 = 1 h 2 I * j-1 f 2 (y)dy, p 1 = 1 h 1 I j f 1 (x)dx, f 2 (y) = y ν-1 f 1 = x λ-ν where h = h 1 = h 2 , I j = [x j-1 , x j ] and I * k = [0, x k ]. Then one computes f 2 -p 2 2 L 2 (0,h) = (ν -1) 2 (1 + (ν -1) 2 )(2(ν -1) + 1) h 2(ν-1)+1 f 2 -p 2 2 L 2 (a,a+h) = a 2(ν-1)+1 η( h a ), a > 0, where η(δ) = (1 + δ) 2(ν-1)+1 -1 2(ν -1) + 1 - [(1 + δ) (ν-1) -1] 2 δ(1 + (ν -1)) 2 and δ > 0. With R * j = ∪ j-1 l=1 R jl = I j × I * j-1 , one notes f 2 -p 2 2 L 2 (I * j-1 ) p 1 2 L 2 (I j ) = h 2(ν-1)+1 (ν -1) 2 (1 + (ν -1) 2 )(2(ν -1) + 1) + j-1 l=1 ( x l h ) 2(ν-1)+1 η( h 2 x l ) x j x j-1 1 h 2 1 x j x j-1 x λ-ν dx 2 dx and x j x j-1 1 h 2 1 x j x j-1 x λ-ν dx 2 dx = x 2λ-2ν+1 j-1 x j-1 h [(1 + h x j-1 ) λ-ν+1 -1] 2 (λ -ν + 1) 2 . As η( h x l ) h 3 x 3 l and 1 δ [(1 + δ) λ-ν+1 -1] 2 (λ -ν + 1) 2 = 1 δ [(λ -ν + 1)δ + (λ -ν + 1)(λ -ν)δ 2 + O(δ 3 )] (λ -ν + 1) 2 = δ + (λ -ν + 1)δ 2 + O(δ 3 ) , we obtain N j=2 f 2 -p 2 2 L 2 (I * j-1 ) p 1 2 L 2 (I j ) = N j=2 h 2(ν-1)+1 [ (ν -1) 2 (1 + (ν -1) 2 )(2(ν -1) + 1) + j-1 l=1 ( x l h ) 2(ν-1)-2 ][(j -1)h] 2λ-2ν+1 (j -1)[(1 + 1 j-1 ) λ-ν+1 -1] 2 (λ -ν + 1)
The estimate

N j=2 f 2 -p 2 2 L 2 (I * j-1 ) p 1 2 L 2 (I j ) = c N j=2 h 2(ν-1)+1 (j -1) 2λ-2ν h 2λ-2ν+1 = h 2(ν-1)+1 h 2λ-2ν+1 N -1 k=1 k 2λ-2ν h 2λ
follows and may be integrated in time. Now consider

f 1 -p 1 2 L 2 (I j ) = x λ-ν -1 h I j x λ-ν dx 2 L 2 (I j ) = a 2(λ-ν)+1 η( h a ), a > 0 (λ-ν) (1+(λ-ν) 2 )(2(λ-ν)+1) h 2(λ-ν)+1 , a = 0 . For j ≥ 2 f 1 -p 1 2 L 2 (I j ) = x 2(λ-ν)+1 j h 2(λ-ν)+1 η h x j h 2(λ-ν)+1 = ( x j h ) 2(λ-ν)-2 h 2(λ-ν)+1 , f 2 2 L 2 (I * j-1 ) = x j-1 0 y 2ν-2 dy = x 2ν-1 j-1 ,
and

f 1 -p 1 2 L 2 (I j ) p 2 2 L 2 (I * j-1 ) = ( x j h ) 2(λ-ν)-2 ( x j-1 h ) 2ν-1 h 2(λ-ν)+1 h 2ν-1 ( x j h ) 2(λ-ν)-2+2ν-1 h 2(λ-ν)+1+2ν-1 . We have N j=2 f 1 -p 1 2 L 2 (I j ) f 2 2 L 2 (I * j-1 ) = c ∞ j=1 j 2(λ-ν)-2+2ν-1 h 2λ-2ν+1+2ν-1 = cζ(3 -2λ)h 2λ .
Again, this estimate may be integrated in time. We now consider case (i): f (t, x, y) = b(t, x)y ν-1 + (χ(x) -1)b(t, x)y ν-1 =:

f 1 + f 2 for ν > 1 2 and Q = [0, 1] 2 = I × I, with I = [0, 1]. Again we define q 1 = Π p t 1 0 b(t, x)dx, q 2 = 1 0 y ν-1 dy . Note that with f 2 (χ(x) -1)b(t, x)y ν-1 ∈ H 0 σ (R + , H 1 (Q)) we get f 2 -q 1 q 2 2 0,0,Q h 2 .
Since

y ν-1 -q 2 2 L 2 (I) ≃ h 2ν-1 , y ν-1 2 L 2 (I) ≤ c and b -q 1 2 0,0,I max{h, ∆t} 2 ∂ t b 2 0,0,I + ∂ x b 2 0,0,I , we have b(t, x)y ν-1 -q 1 (t, x)q 2 (y) 2 0,0,Q b 2 0,0,I , y ν-1 -q 2 2 L 2 (I) + y ν-1 2 L 2 (I) , b -q 1 2 0,0,I (26)
h 2ν-1 + max{h, ∆t} 2 .
For the lower bound, we first consider the error in R 11 resulting from the corner singularity. There the error of approximation by a spatially constant function c is given by

a(t, θ)r λ-1 -c(t) 2 0,0,R 11 ∞ 0 h 0 π/2 0 r(a(t, θ)r λ-1 -c) 2 dθ drd σ t .
If a = 0, we may find a small intervall (θ 0δ, θ 0 + δ) on which a is nonzero for a time interval I. We estimate

a(t, θ)r λ-1 -c 2 0,0,R 11 I h 0 θ 0 +δ θ 0 -δ r(a(t, θ)r λ-1 -c) 2 dθ dr d σ t .
As a is a restriction of an eigenfunction of the Laplace-Beltrami operator, it is smooth, and up to higher order terms (h.o.t.) in h we compute

I h 0 θ 0 +δ θ 0 -δ r(a(t, θ)r λ-1 -c) 2 dθ dr = h 0 θ 0 +δ θ 0 -δ r(a(t, θ 0 )r λ-1 -c) 2 dθ dr d σ t + h.o.t. = 2δ I h 0 r(a(θ 0 )r λ-1 -c) 2 dr d σ t + h.o.t.
Now we may explicitly compute the infimum over c:

I h 0 r(a(t, θ 0 )r λ-1 -c) 2 dr ≥ C σ |I| min t∈I a(t, θ 0 ) 2 (λ -1) 2 2λ(λ + 1) 2 h 2λ .
This lower order bound for the error of order h λ matches with the upper bound from above.

An analogous argument for the regular edge intensity factor, case (i), shows that the inequality ( 26) is an equality up to higher order terms in h. We therefore obtain b(x)y ν-1q 1 (x)q 2 (y) 0,0,Q h ν-1 2 + h.o.t.

In a neighborhood of the edge or corner, the solution is given by its singular expansion, and the approximation error coincides with the approximation error for the singular functions,

Here, for every triangle △ and every time interval I n = [t n-1 , t n ] we define the partial error indicators

η △,∇ Γ (I n ) 2 = h △ tn t n-1 △ [∇ Γ ( ḟ -V φh,∆t )] 2 ds x dt , η △,∂t (I n ) 2 = ∆t tn t n-1 △ [∂ t ( ḟ -V φh,∆t )] 2 ds x dt .
The time integral is approximated by the trapezoidal rule, and the tangential gradient of a function F is computed as

∇ Γ F (t, x) = P Γ ∇F = ∇F (t, x) -ν(ν • ∇F (t, x))
with the outer unit normal vector ν to Γ, resp. the projection P Γ onto the tangent bundle of Γ.

To compute η △,∇ Γ from φh,∆t , we consider the gradient of V φh,∆t as a singular integral:

∇ Γ V φh,∆t (t, x) = -1 4π P Γ Γ (x -y) φh,∆t (t -|x -y|, y) |x -y| 3 + φh,∆t (t -|x -y|, y) |x -y| 2 ds y = -1 4π Nt m=1 Ns i=1 ϕ m i P Γ Γ ϕ i (y) β m (t -|x -y|) x -y |x -y| 3 + βm (t -|x -y|)
xy |x -y| 2 ds y .

Using the explicit form of β m , we obtain

∇ Γ V φh,∆t (t, x) = -1 4π Nt m=1 Ns i=1 ϕ m i t -t m-1 △t P Γ t-tm≤|x-y|≤t-t m-1 ξ i (y) x -y |x -y| 3 ds y - t -t m+1 △t P Γ t-t m+1 ≤|x-y|≤t-tm ξ i (y) x -y |x -y| 3 ds y .
The integrals are evaluated with a composite hp-graded quadrature, like the entries of the BEM matrix in [START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF]. See [START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF] for details.

Numerical experiments

Example 1: We consider the Dirichlet problem Vφ = f on the unit sphere Γ = S 2 with the right hand side f (t, x, y, z) = sin(t) 5 x 2 and [0, T ] = [0, 2.5]. We use a discretization by linear ansatz and test functions in space and time. Γ is approximated by uniform meshes of 80, 320,1280, and 5120 triangles, and the time step ∆t is 0.4, 0.2, 0.1, resp. 0.05 for the respective meshes to keep ∆t h fixed. The numerical results are compared to the exact solution.

Figure 1 shows the convergence of the error φφ h,∆t in the energy norm as well as the L 2 error in the sound pressure and compares them to both residual and ZZ error indicators. The resulting convergence rates are similar: We obtain a convergence rate of 0.93 in energy norm, 0.97 in sound pressure, 0.9 in the residual error indicator, and 1.02 in the ZZ indicator. This illustrates the reliability and efficiency of both error indicators with respect to the energy norm and related quantities such as the sound pressure in an example with known exact solution. More precisely, the quotient of the error estimate and the energy error, the efficiency index, remains approximately constant at 0.025 as the number of degrees of freedom increases.

Example 2:

We consider the Dirichlet problem Vφ = f on the square screen Γ = [-0.5, 0.5] 2 × {0} with the right hand side f (t, x, y, z) = sin(t) 5 x 2 for times [0, 2.5]. Using a discretization by linear ansatz and test functions in space and time, we compare the error of a uniform discretization to the error of an adaptive series of meshes, steered by the residual error estimate. The time step is fixed at ∆t = 0.1, and the uniform meshes consist of 18, 288, 648, 1352, and 6050 triangles, while the adaptive refinements correspond to 36, 74, 164, 370, 784, 1676, 3485, and 7432 triangles.

Figure 2 shows the convergence of the error indicator and the error in the energy norm, for both the uniform and adaptive series of meshes. The convergence rate is approximately 0.48 for uniform refinements, compared to 0.77 for adaptive refinements. The convergence rate in the uniform case agrees with the theoretical prediction of 0.5 from [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF], and the adaptive convergence rate of 0.77 recovers the results for time-independent screen problems [START_REF] Carstensen | A posteriori error estimate and h-adaptive algorithm on surfaces for Symm's integral equation[END_REF].

As in the elliptic case, the convergence rate of the adaptive refinements does not reach the optimal rate of 1.5 achieved with algebraically graded meshes, as demonstrated in [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF]. The optimal anisotropic graded meshes cannot be obtained by mesh refinements: While adaptive meshes are locally quasi-uniform, graded meshes involve arbitrarily thin triangles with shallow angles near the edges of the screen. A heuristic explanation for the substantially higher rates of (anisotropic) graded meshes is contained in [START_REF] Carstensen | Residual-based a posteriori error estimate for hypersingular equation on surfaces[END_REF]. Figure 4 shows the convergence of the error indicator and the error in the energy norm, for both the uniform and adaptive series of meshes. The convergence rate is approximately 0.49 for uniform refinements, compared to 0.78 for adaptive refinements, almost identical to the square screen in Example 2.

Figure 5 shows representative adaptive meshes, where the color scale highlights the residual-based indicator values for each element. As expected, mesh refinements concentrate in the two sharper corners of the triangle.

Example 4: We consider the Dirichlet problem Vφ = f on the triangle Γ with angles of 30, 60 and 90 degrees, as depicted in Figure 7. The right hand side is given by f (t, x, y, z) = sin(t) 5 , and we consider times [0, 2.5]. Using the discretization from Example 2, we compare the error on uniform meshes to the error of an adaptive series of meshes, steered by the residual error estimate. The time step is fixed at ∆t = 0.1.

Figure 6 shows the convergence of the error indicator and the error in the energy norm, for both the uniform and adaptive series of meshes. The convergence rate is approximately 0.448 for uniform refinements, compared to 0.65 for adaptive refinements. The rates are slightly reduced compared to Examples 2 and 3, possibly because the asymptotic regime only sets in for higher degrees of freedom because of the small angles of 30 degrees in the triangulation. From Experiments 2, 3 and 4 we conclude that the convergence rate is does not depend on the angles of the triangle, and therefore the corner singularity. The convergence rate of around 1 2 on uniform meshes matches the rate theoretically expected for the approximation of the edge singularity [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF], while the approximation error from the corner singularities is of higher order. The adaptive convergence rates of around 0.78 are compatible with the convergence rates of around 0.8 for the time-independent Laplace equation in [START_REF] Carstensen | A posteriori error estimate and h-adaptive algorithm on surfaces for Symm's integral equation[END_REF]. The rates are slightly reduced in Example 4, with angles of 30 degrees, possibly because of the

2 ,u

 2 and ||u|| r,s,σ,Γ, * < ∞} . D ′ + (E) respectively S ′ + (E) denote the spaces of distributions, respectively tempered distributions, on R with support in [0, ∞), taking values in E = H s (Γ), H s (Γ). The relevant norms are given by u r,s,σ := u r,s,σ,Γ = +∞+iσ -∞+iσ |ω| 2r û(ω) 2 s,ω,Γ dω 1 r,s,σ, * := u r,s,σ,Γ, * = +∞+iσ -∞+iσ |ω| 2r û(ω) 2 s,ω,Γ, * dω 1 2

Figure 1 :

 1 Figure 1: Energy error, residual and ZZ error indicators for Dirichlet problem on Γ = S 2 , Example 1.

Figure 2 :

 2 Figure 2: Energy error and residual error indicators for Dirichlet problem on Γ = [-0.5, 0.5] 2 × {0}, Example 2.

Figure 3

 3 Figure3shows representative adaptive meshes, where the color scale highlights the residual-based indicator values for each element. Mesh refinements concentrate at the left and right edges, where the right hand side is steep, and to a lesser extent also at the top and bottom edges.

Example 3 :

 3 We consider the Dirichlet problem Vφ = f on the triangle Γ with angles of 45, 45 and 90 degrees, as depicted in Figure5. The right hand side is given by f (t, x, y, z) = sin(t)5 , and we consider times [0, 2.5]. Using the discretization from Example 2, we compare the error on uniform meshes to the error of an adaptive series of meshes, steered by the residual error estimate. The time step is fixed at ∆t = 0.1.

Figure 3 :

 3 Figure 3: Meshes 1, 2, 3 and 6 generated by adaptive refinements, Example 2.

Figure 4 :

 4 Figure 4: Energy error and residual error indicators for Dirichlet problem isosceles triangle, Example 3.

Figure 5 :

 5 Figure 5: Meshes 3, 5, 7 and 8 generated by adaptive refinements, Example 3.

Figure 7

 7 Figure 7 shows representative adaptive meshes, where the color scale highlights the residual-based indicator values for each element. As expected, mesh refinements concentrate in the corners according to their sharpness.

Figure 6 :

 6 Figure 6: Energy error and residual error indicators for Dirichlet problem on 30-60-90 triangle, Example 4.

Figure 7 :

 7 Figure 7: Meshes 2, 5, 8 generated by adaptive refinements, Example 4.

  with the same norm to a bounded operator on H -s ω (Γ). If the triangulation is globally quasi-uniform, we show that Q h is a bounded operator on H 1

ω (Γ) with operator norm uniformly bounded in ω: Lemma 2.5. Assume that Γ does not have reentrant corners. Consider a globally quasiuniform mesh and s ∈ [0, 1]. For all φ ∈ H s ω (Γ), we have:

H. G. acknowledges partial support by ERC Advanced Grant HARG 268105 and the EPSRC Impact Acceleration Account. C. Ö. is supported by a scholarship of the Avicenna Foundation.

b(t, x)y ν-1 , y λ-ν i ρ ν-1 , respectively a(t, θ)r λ-1 , up to lower order terms in h. We conclude that E(φ, h, ∆) = φ -Π h,∆t φ 2,0,Γ max{h, ∆t} min{ν-1 2 ,λ} + h.o.t. .

References [START_REF] Petersdorff | Randwertprobleme der Elastizitätstheorie für Polyeder-Singularitäten und Approximation mit Randelementmethoden[END_REF][START_REF] Petersdorff | Regularity of mixed boundary value problems in R 3 and boundary element methods on graded meshes[END_REF] show how to deduce approximation results for piecewise linear functions on triangular meshes from piecewise bilinear functions on rectangles. Altogether, the proof of Theorem 6.1 is complete.

Algorithmic details

The a posteriori error estimate from Theorem A leads to an adaptive mesh refinement procedure, based on the four steps:

The precise algorithm is given as follows:

Adaptive Algorithm: Input: Spatial mesh T = T 0 , refinement parameter θ ∈ (0, 1), tolerance ǫ > 0, data f .

1. Solve V φh,∆t = ḟ on T .

2. Compute the error indicators η(△) in each triangle △ ∈ T .

3. Find η max = max △ η(△).

Stop if

5. Mark all △ ∈ T with η(△ i ) > θη max .

6. Refine each marked triangle into 4 new triangles to obtain a new mesh T Choose ∆t such that ∆t ∆x ≤ 1 for all triangles. 7. Go to 1.

Output: Approximation of φ.

In the first step, we solve V φ = ḟ using the Galerkin discretization [START_REF] Gimperlein | Adaptive time domain boundary element methods and engineering applications[END_REF]

where β m is the piecewise linear hat function in time associated to time t m ,

and ξ i is the piecewise linear hat function in space associated to node i. As a step towards adaptive mesh refinements in space-time, we here focus on timeintegrated error indicators as they are relevant for geometric singularities. As shown in [START_REF] Gimperlein | Boundary elements with mesh refinements for the wave equation[END_REF], for polyhedral meshes and screens time-independent graded meshes lead to quasioptimal convergence rates in spite of singularities of the solutions. The time-integrated error indicator for triangle △ is computed as


Tibor Illés

Marianna Nagy

Tamás Terlaky

Polynomial Interior Point Algorithms for General Linear Complementarity Problems

Keywords: linear complementarity problem, sufficient matrix, P * -matrix, interior point method, affine scaling method, predictor-corrector algorithm

Linear Complementarity Problems (LCPs) belong to the class of NP-complete problems. Therefore we can not expect a polynomial time solution method for LCPs without requiring some special property of the coefficient matrix. Following our recently published ideas we generalize affine scaling and predictor-corrector interior point algorithms to solve LCPs with general matrices in EP-sense, namely, our generalized interior point algorithms either solve the problems with rational coefficient matrix in polynomial time or give a polynomial size certificate that our matrix does not belong to the set of P * (κ) matrices, with arbitrary large, but apriori fixed, rational, positive κ.

Introduction

Consider the Linear Complementarity Problem (LCP): find vectors x, s ∈ R n that satisfy -M x + s = q, x s = 0, x, s ≥ 0, [START_REF] Dikin | Iterative solution of problems of linear and quadratic programming[END_REF] where M ∈ R n×n and q ∈ R n , and the notation xs is used for the coordinatewise (Hadamard) product of the vectors x and s.

LCPs belong to the class of NP hard problems, since the feasibility problem of linear equations with binary variables can be described as an LCP [START_REF] Kojima | A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems[END_REF]. Therefore we can not expect an efficient (polynomial time) solution method for LCPs without requiring some special property of the matrix M .

⋆ The research of Tibor Illés and Marianna Nagy has been supported by the Hungarian National Research Fund OTKA No. T049789. The work of T. Terlaky has been supported by a grant from Lehigh University. Email: Tibor Illés [tibor.illes@strath.ac.uk], Marianna Nagy [nmariann@cs.elte.hu], Tamás Terlaky [terlaky@lehigh.edu].

In [START_REF] Illés | A polynomial path-following interior point algorithms for general linear complementarity problems[END_REF] we modified long-step path-following Interior Point Methods (IPMs) for LCPs with a general coeffi- cient matrix M . The modified algorithm either solves the LCP, or gives a certificate that the matrix of the problem is not P * (κ) (with apriori given but arbitrary large κ), or gives a certificate that the LCP has no solution. Algorithms that do not require any special property of the Matrix M are needed becasue it cannot be verified in polynomial time if matrix M belongs to the class of matrices that allow polynomial time solvability of the LCP. Indeed, Tseng [START_REF] Tseng | Co-NP-completeness of some matrix classification problems[END_REF], proved that the problem deciding whether a square matrix with rational entries is a column sufficient matrix is co-NP-complete, suggesting that it can not be decided in polynomial time whether there is a finite nonnegative κ with which matrix M is P * (κ). In this paper we show, that the idea behind the modification of long-step path-following algorithm is general, i.e., it can be adapted to other IPMs as well. Thus here we present the modifications of two well known IPMs: the affine scaling and predictor corrector algorithms.

In [START_REF] Illés | A polynomial path-following interior point algorithms for general linear complementarity problems[END_REF] first we discussed the generalization of the embedding technique of Kojima et al. [START_REF] Kojima | A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems[END_REF]. The embedding technique is one of the options to ensure the availability of an initial interior point. This technique is independent of the particular IPM, therefore we do not repeat it in this paper. Thus we may assume without loss of generality that an initial interior feasible solution is known.

The rest of the paper is organized as follows. The next section deals with the fundamental properties of P * (κ)-matrices and some well-known results are presented. In Section 3 we summarize the results of paper [START_REF] Illés | A polynomial path-following interior point algorithms for general linear complementarity problems[END_REF]. Section 4 deals with the modification of two well known interior point algorithms: the affine scaling and predictor-corrector algorithms.

For ease of understanding and for self containedness the main results of the papers [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF][START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF] are summarized in the Appendix.

Notation:

We use the following notations throughout the paper. Scalars and indices are denoted by lowercase Latin letters, vectors by lowercase boldface Latin letters, matrices by capital Latin letters, and finally sets by capital calligraphic letters. Let R n ⊕ (R n + ) denote the nonnegative (positive) orthant of R n . Further, I denotes the identity matrix of appropriate dimension, and X is the diagonal matrix whose diagonal elements are the coordinates of the vector x, so X = diag(x). The vector x s = Xs is the componentwise product (Hadamard product) of the vectors x and s, and for α ∈ R the vector x α denotes the vector whose ith component is x α i . The largest and smallest coordinate of a vector is denoted by max(x) and min(x), respectively. We denote the vector of ones by e. Furthermore, for a vector x and a matrix M we define the sets

I + (x) = {1 ≤ i ≤ n : x i (M x) i > 0} and I -(x) = {1 ≤ i ≤ n : x i (M x) i < 0}, which are used in the definition of P * (κ) matrices. Finally, F 0 := (x, s) ∈ R 2n
+ : -M x + s = q denotes the set of strictly feasible solutions of the LCP.

Matrix classes and the Newton step

The class of P * (κ)-matrices were introduced by Kojima et al. [START_REF] Kojima | A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems[END_REF], and it can be considered as a generalization of the class of positive semidefinite matrices.

Definition 1 Let κ ≥ 0 be a nonnegative number. A matrix M ∈ R n×n is a P * (κ)-matrix if for all x ∈ R n (1 + 4κ) i∈I + (x) x i (M x) i + i∈I -(x) x i (M x) i ≥ 0.(2)
The nonnegative real number κ denotes the weight that need to be used at the positive terms so that the weighted 'scalar product' is nonnegative for each vector x ∈ R n . Therefore, naturally P * (0) is the class of positive semidefinite matrices (setting aside the symmetry of the matrix M ). Definition 2 A matrix M ∈ R n×n is called a P *matrix if it is a P * (κ)-matrix for some κ ≥ 0, i.e.

P * = κ≥0 P * (κ).
The class of sufficient matrices was introduced by Cottle, Pang and Venkateswaran [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF].

Definition 3 A matrix M ∈ R n×n is a column suffi- cient matrix if for all x ∈ R n X(M x) ≤ 0 implies X(M x) = 0,
and row sufficient if M T is column sufficient. Matrix M is sufficient if it is both row and column sufficient.
Kojima et al. [START_REF] Kojima | A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems[END_REF] proved that any P * -matrix is column sufficient and Guu and Cottle [START_REF] Guu | On a subclass of P0[END_REF] proved that it is row sufficient, too. Therefore, each P * -matrix is sufficient. Väliaho proved the other direction of inclusion [START_REF] Väliaho | P * -matrices are just sufficient[END_REF], thus the class of P * -matrices coincides with the class of sufficient matrices. Definition 4 A matrix M ∈ R n×n is a P 0 -matrix, if all of its principal minors are nonnegative.

For further use we recall some results about P * (κ)and P 0 -matrices. The reader may consult the book of Kojima et al. [START_REF] Kojima | A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems[END_REF]Lemma 4.1 p. 35] for the proof of the following proposition. Proposition 5 A matrix M ∈ R n×n is a P 0 -matrix if and only if

M ′ = -M I S X
is a nonsingular matrix for any positive diagonal matrices X, S ∈ R n×n . 2 Proposition 5 enables us to check in strongly polynomial time whether matrix M is P 0 or not. The next statement is used to guarantee the existence and uniqueness of Newton directions that are the solution of system (3) for various values of vector a ∈ R n , where a depends on the particular IPM. Corollary 6 Let M ∈ R n×n be a P 0 -matrix, (x, s) ∈ F 0 . Then, for all a ∈ R n the system

-M ∆x + ∆s = 0 s∆x + x∆s = a (3)
has a unique solution (∆x, ∆s). 2

The following estimations for the Newton direction are used in the complexity analysis of IPMs. The next lemmas are proved by Potra in [START_REF] Potra | The Mizuno-Todd-Ye algorithm in a larger neighborhood of the central path[END_REF]. Lemma 7 Let (x, s) ∈ F 0 and M be an arbitrary n×n real matrix and (∆x, ∆s) be a solution of system [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF]. Then

i∈I + ∆x i ∆s i ≤ 1 4 a √ xs 2 .
Lemma 8 Let the matrix M be a P * (κ)-matrix, x, s ∈ F 0 , and a ∈ R n . Let (∆x, ∆s) be the solution of system [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF]. Then

∆x∆s ∞ ≤ 1 4 + κ a √ xs 2 , ∆x∆s 1 ≤ 1 2 + κ a √ xs 2 , ∆x∆s 2 ≤ 1 4 + κ 1 2 + κ a √ xs 2 .
The first statement's proof in the previous lemma is similar to the proof of Lemma 5.1 by Illés, Roos and Terlaky [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF]. The second estimation follows from the Lemma 7 by using some properties of P * (κ)-matrices, and the last estimation is a corollary of the first and second statements using some properties of norms.

Let the current point be (x, s) ∈ F 0 and (∆x, ∆s) be the corresponding Newton direction. The new point with step length θ is given by (x(θ), s(θ)) = (x + θ∆x, s + θ∆s). We use the following notations for scaling

v = xs µ , g = ∆x∆s µ , (4)
where in the affine scaling algorithm for the purpose of scaling we have µ ≡ 1, otherwise µ > 0. In affine scaling algorithms we use the δ a centrality measure

δ a (xs) = max √ xs min √ xs .
In the predictor-corrector algorithm the so-called negative infinity neighborhood D(γ), defined by Potra and Liu in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF], is used. The negative infinity neighborhood for γ ∈ (0, 1) is defined as

D(γ) := (x, s) ∈ F 0 : x s ≥ γ x T s n ,
where F 0 := (x, s) ∈ R 2n + : -M x + s = q denotes the set of strictly feasible solutions of the LCP. The D(γ) neighborhood is considered to be a "wide neighborhood".

Previous results

In this section we restate those results of paper [START_REF] Illés | A polynomial path-following interior point algorithms for general linear complementarity problems[END_REF], that we use in the rest of this paper.

The inequality in the definition of P * (κ)-matrices gives the following lower bound on κ for any vector

x ∈ R n : κ ≥ κ(x) = - 1 4 x T M x i∈I + x i (M x) i .
The following two lemmas are immediate consequences of the definition of P * (κ) and P * -matrices. Lemma 9 Let M be a real n × n matrix and κ > 0 be a given parameter. If there exists a vector x ∈ R n such that κ(x) > κ, then the matrix M is not P * (κ) and x is a certificate for this fact. Lemma 10 Let M be a real n × n matrix. If there exists a vector x ∈ R n such that I + (x) = {i ∈ I :

x i (M x) i > 0} = ∅, then the matrix M is not P * and
x is a certificate for this fact.

Our modified interior point algorithms are based on the assumption that the problem data is rational, and we are interested either to provide a polynomial size solution to the problem, or to provide a certificate that our matrix M is not a P * (κ)-matrix for some κ > 0. During all iterations of our algorithm we have an actual estimate for κ, say κ. At the beginning of our algorithms, we assume that κ = 0. In every iteration at different stages of the algorithm we are checking whether our assumption related to κ is violated. If we detect, that the matrix M is not P * (κ) with actual κ, then the value of κ will be increased. Its new value will be the lower bound defined by the actual Newton direction ∆x. In IPMs the P * (κ) property need to hold only for the actual Newton direction ∆x in various ways, for example this property ensures that with a certain step size the new iterate is in an appropriate neighborhood of the central path and/or the complementarity gap is sufficiently reduced. Consequently, if the desired results do not hold with the current κ value, we update κ by increasing it to the lower bound determined by the Newton direction ∆x, i.e.,

κ = κ(∆x) =: -∆x T ∆s 4 i∈I + ∆x i ∆s i (∆s = M ∆x). (5)
The following lemma is our main tool to verify when matrix M fails to satisfy the P * (κ) property. Furthermore, the concerned vector ∆x is a certificate, whose encoding size is polynomial when it is computed as the solution of the Newton system (3) from rational data.

We use this lemma during the analysis. The first statement is simply the negation of the definition. We point out in Lemma 12 that if Lemma 4.3 of [START_REF] Jansen | A family of polynomial affine scaling algorithms for positive semidefinite linear complementarity problems[END_REF] does not hold, then the second statement is realized. We show in Lemma 15 and Lemma 16 that if Theorem 3.3 of [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF] does not hold then the second, the third or the fourth statement is realized. Lemma 11 Let M be a real n × n matrix, κ ≥ 0 be a given parameter, and (x, s) ∈ F 0 . If any of the following statements holds then the matrix M is not a P * (κ)-matrix.

(1) There exists a vector y ∈ R n such that

(1 + 4κ) i∈I + (y) y i w i + i∈I -(y) y i w i < 0,
where w = M y and I + (y) = {i ∈ I :

y i w i > 0}, I -(y) = {i ∈ I : y i w i < 0}.
(2) There exists a solution (∆x, ∆s) of system [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF] such that

∆x∆s ∞ > 1 + 4κ 4 a √ xs 2 .
(3) There exists a solution (∆x, ∆s) of system [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF] such that

max   i∈I + ∆x i ∆s i , - i∈I - ∆x i ∆s i   > 1 + 4κ 4 a √ xs 2 .
(4) There exists a solution (∆x, ∆s) of system [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF] such that

∆x T ∆s < -κ a √ xs 2 .

Interior point algorithms in EP form

Hereafter we modify the following two popular families of IPMs:

• A family of affine scaling algorithms [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF]:

The Newton direction is the solution of system [START_REF] Cottle | Sufficient matrices and the linear complementarity problem[END_REF] with µ = 1 and a = -v 2r+2 v 2r , where r ≥ 0 is the degree of the algorithm.

The choice of a implies a √ xs

2 = v 2r+1 2 v 2r 2 .
• Predictor-corrector algorithms [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]:

The predictor Newton direction is the solution of system (3) with a = -xs (the affine scaling direction with r = 0).

The choice of a implies a √ xs

2 = x T s;
The corrector (centering) Newton direction is the solution of system (3) with a = µexs, where

µ = x T s n . The choice of a implies a √ xs 2 = µ xs µ - µ xs 2 .

Affine scaling IPMs

First we deal with affine scaling IPMs. We modify the family of algorithms proposed in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF], where the particular algorithms correspond to the degree r ≥ 0 of the algorithm, where r = 0 gives the classical primaldual affine scaling algorithm, while r = 1 gives the primal-dual Dikin affine scaling algorithm [START_REF] Dikin | Iterative solution of problems of linear and quadratic programming[END_REF]. Further, there is a step length parameter ν, that depends on the degree r (defined among the inputs of the algorithm), and µ ≡ 1 in scaling (4).

We check not only the solvability and uniqueness of the Newton system, but also the decrease of the complementarity gap after a step. For the actual value of κ we determine θ * a (κ), which is a theoretical lower bound for the maximal feasible step length in the specified neighborhood if the matrix M satisfies the P * (κ) property. Therefore, if after a step the decrease of the complementarity gap is not large enough, it means, that the matrix M is not P * (κ) with the actual value of κ, so we update κ or exit the algorithm with a proper certificate. If the new value of κ can not be defined by ( 5), then the matrix M is not P * , so we stop and the Newton direction ∆x is a certificate. If the new value of κ is larger than κ, then the matrix is not P * (κ), therefore the algorithm stops as well and ∆x is a certificate. In the rest of this subsection we consider the case r > 0. The modified algorithm is presented in Algorithm 1.

Illés et al. proved [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF], that if the matrix M is a P * (κ)-matrix, then the step length θ * a (κ) is feasible, with that step size the new iterate stays within the specified neighborhood and it provides the required decrease of the complementarity gap. The following lemma shows, if the decrease of the complementarity

Algorithm 1 -Affine scaling algorithm

Input:

an upper bound κ > 0 on the value of κ; an accuracy parameter ε > 0; a centrality parameter τ ; the degree of scaling r > 0; a strictly feasible initial point (x 0 , s 0 ) ∈ F 0 such that δ a (x 0 s 0 ) ≤ τ ;

ν := 2/ √ n, if 0 < r ≤ 1 2τ 2-2r / √ n, if 1 ≤ r; θ * a (κ) := min 2 (1 + 4κ)τ 1 + 4κ + 1 τ 2 n - 1 τ √ n , √ n (r + 1)τ 2r , 4(τ 2r -1) (1 + 4κ)(1 + τ 2 )τ 2r √ n , ν .

begin

x := x 0 , s := s 0 , κ := 0; while x T s ≥ ε do calculate the Newton direction (∆x, ∆s) with a = -v 2r+2 / v 2r ; if (the Newton direction does not exist or it is not unique) then return the matrix is not P 0 ;

% see Corollary 6 θ = argmin x(θ) T s(θ) : δ a (x(θ), s(θ)) ≤ τ, (x(θ), s(θ)) ≥ 0 ; if x( θ) T s( θ) > (1 -0.25 ν θ * a (κ))
x T s then calculate κ(∆x); % see [START_REF] Fukuda | EP theorems and linear complementarity problems[END_REF] if (κ(∆x) is not defined) then return the matrix is not P * ;

% see Lemma 10 if (κ(∆x) > κ) then return the matrix is not P * (κ);

% see Lemma 9 κ = κ(∆x); update θ * a (κ);

% it depends on κ

x = x( θ), s = s( θ); end end. gap is not sufficient, then the matrix M does not belong to the class of P * (κ)-matrices.

Lemma 12 If x( θ) T s( θ) > (1 -0.25 ν θ * a (κ))
x T s, that is, the decrease of the complementarity gap within the δ a ≤ τ neighborhood is not sufficient, then the matrix M of the LCP is not P * (κ) with the actual value of κ. The Newton direction ∆x serves as a certificate.

Proof: Based on [12, Lem 4.3] (see Lemma 20 in the Appendix) the complementarity gap at θ * a (κ) is smaller than (1 -0.25 ν θ * a (κ)) x T s, furthermore by Theorem 21, if M is a P * (κ)-matrix, then the point

(x * , s * ) = (x(θ * a (κ)), s(θ * a (κ))) is feasible. Therefore, if x( θ) T s( θ) > (1 -0.25 ν θ * a (κ))
x T s, then because the step length θ * a (κ) is not considered in definition of θ (see the affine scaling algorithm), so either (x * , s * ) is not feasible, or this point is not in the τ neighborhood of the central path, namely δ a (x * s * ) > τ . We show that both cases imply, that the matrix M is not P * (κ) with the actual κ value.

Let us denote the first three terms in the definition of θ * a (κ) by θ 1 , θ 2 , and θ 3 , respectively. We follow the proof of Theorem 6.1 in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF] (see Theorem 21 in the Appendix). We need to reconsider only the expressions depending on κ. There exist positive constants α and β such that β α = τ 2 and αe ≤ v 2 ≤ βe. The function

ϕ(t) = t -θ t r+1 v 2r remains monotonically increasing on the interval [0, β] if θ ≤ θ 2 ≤ v 2r
(r+1)β r . Additionally, inequalities [START_REF] Potra | The Mizuno-Todd-Ye algorithm in a larger neighborhood of the central path[END_REF] in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF] hold for θ ≤ θ 2 , thus for θ * a (κ) too:

min(v * 2 ) ≥ α -θ * a (κ) α r+1 v 2r -(θ * a (κ)) 2 g ∞ , (6)
max(v * 2 ) ≤ β -θ * a (κ) β r+1 v 2r + (θ * a (κ)) 2 g ∞ , ( 7
)
where g is defined by (4) (see p.3). Let us first consider the case δ a (x * s * ) > τ , i.e., max(x * s * ) > τ 2 min(x * s * ). From inequalities ( 6) and ( 7) one has

τ 2 α -θ * a (κ) α r+1 v 2r -(θ * a (κ)) 2 g ∞ < β -θ * a (κ) β r+1 v 2r + (θ * a (κ)) 2 g ∞ .
Recalling the equality ατ 2 = β, and dividing both sides of the inequality by θ * a (κ), one gets

β r -α r v 2r < θ * a (κ) 1 α + 1 β g ∞ . (8)
If θ * a (κ) is substituted by θ 3 = 4(τ 2r -1) (1+4κ)(1+τ 2 )τ 2r √ n , the right hand side of inequality ( 8) increases, so the inequality is still true. After substitution one has

1 + 4κ 4 β < g ∞ . (9)
Since v 2 ≤ βe and g = ∆x∆s (see the notation given by ( 4)), inequality (9) gives

v 2 ∞ ≤ β < 4 1 + 4κ g ∞ = 4 1 + 4κ ∆x∆s ∞ . ( 10
)
One can check, that

v 2r+1 2 ≤ v 2 ∞ v 2r 2 . ( 11
)
Since (∆x, ∆s) is the solution of system (3) with a = -v 2r+2 / v 2r , by inequalities [START_REF] Illés | An EP theorem for dual linear complementarity problems[END_REF] and [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF] we have

a √ xs 2 = v 2r+1 v 2r 2 ≤ v 2 ∞ .
Therefore, by the second statement of Lemma 11, we get that inequality (10) contradicts to the P * (κ) property and vector ∆x is a certificate for this fact. Now we consider the case (x * , s * ) is not feasible, so there exists such an index i, that either x * i < 0 or s * i < 0. Let us consider the maximum feasible step size θ < θ * a (κ), for which (x( θ), s( θ)) = (x, ŝ) ≥ 0 holds and at least one of its coordinates is 0. For this point x ŝ = 0, else θ = θ by the definition of θ, and the new point would be an exact solution, so the decrease of the complementarity gap would be x T s contradicting with the assumption of the lemma. Therefore 0 = max(x ŝ) > τ 2 min(x ŝ) = 0, so inequality (8) holds with θ. Because of θ 3 ≥ θ * a (κ) > θ, inequality (9) holds as well, and as we have already seen this means that the matrix M is not P * (κ) and the vector ∆x is a certificate for this fact.

The following lemma proves, that the algorithm is well defined. Lemma 13 At each iteration, when the value of κ is updated, then the new value of θ * a (κ) satisfies the inequality x( θ) T s( θ) ≤ (1 -0.25 ν θ * a (κ))x T s. Proof: In the proof of Theorem 6.1 in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF] (see Theorem 21 in the Appendix) the P * (κ) property is only used for the vector ∆x. When parameter κ is updated, then we choose the new value in such a way, that the inequality in the definition of P * (κ)-matrices (2) would hold for vector ∆x. Therefore the new point defined by the updated value of step size θ * a (κ) is feasible and it is in the τ -neighborhood of the central path. Thus the new value of θ * a (κ) was considered in the definition of θ, so x( θ) T s( θ) ≤ (1 -0.25 ν θ * a (κ)) x T s. Now we are ready to state the complexity result for the modified affine scaling algorithm for general LCPs in case an initial interior point is given.

Theorem 14 Let (x 0 , s 0 ) ∈ F 0 such that δ a (x 0 s 0 ) ≤ τ = √ 2.
Then after at most p iterations the affine scaling algorithm either yields a vector (x, ŝ) such that xT ŝ ≤ ε and δ a (xŝ) ≤ τ , or it gives a polynomial size certificate that the matrix is not P * (κ), where κ ≤ κ is the largest value of parameter κ and

p =                O n(1+4κ) 1-2 -r log (x 0 ) T s 0 ε , if 0 < r ≤ 1, n ≥ 4 O n(1 + 4κ) log (x 0 ) T s 0 ε , if r = 1, n ≥ 4 O 2 2r-2 n(1 + 4κ) log (x 0 ) T s 0 ε , if 1 <
r and n sufficiently large.

Proof:

The algorithm at each iteration either takes a step, or detects, that the matrix is not P * (κ) and stops. If we take a Newton step, then by the definition of the algorithm and by Lemma 13 the decrease of the complementarity gap is at least 0.25 ν θ * a (κ) x T s. One can see from the definition of θ * a (κ) that larger κ means smaller θ * a (κ), so smaller lower bound on the decrease of the complementarity gap. Therefore, if the algorithm stops with an ε-optimal solution, then each Newton step decreases the complementarity gap by more than 0.25 ν θ * a (κ) x T s. It means that after at most as many steps as in the original method the complementarity gap decreases below ε in case for each vector during the algorithm sufficient decrease of the complementarity gap is realized according to the P * (κ) property or at an earlier iteration the lack of P * (κ)-property is detected. This observation, combined with the complexity theorem of the original algorithm, see Theorem 6.1 in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF] (see also Theorem 22 in the Appendix) proves our statement.

At the end of this subsection let us note that the case r = 0 can be treated analogously.

Predictor-corrector IPMs

In this section we modify the algorithm proposed in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]. In this predictor-corrector algorithm we take affine and centering steps alternately. In a predictor step θ * p (κ) (see the definition in Lemma 15) is a theoretically feasible step length if the matrix M is P * (κ). Therefore, if the maximal feasible step length is smaller than θ * p (κ), then the matrix is not P * (κ) with the actual value of κ, so κ should be increased. In a corrector step we return to the smaller D(γ) neighborhood with step size θ * c (κ) (see the definition in Lemma 16) if the matrix is P * (κ). Accordingly, if the new point with step length θ * c (κ) is not in D(γ), then the matrix M is not P * (κ) with the actual value of κ, so κ should be updated. Similarly to the Affine Scaling algorithm, if in a predictor or corrector step the new value of κ cannot be defined by [START_REF] Fukuda | EP theorems and linear complementarity problems[END_REF], then the matrix is not P * and the current Newton direction is a certificate of it. Furthermore, if the new value of κ is larger than κ, then the matrix is not P * (κ) and the Newton direction is a certificate for it. The modified algorithm is summarized in Algorithm 2.

Potra and Liu determined the maximum feasible predictor step length as the minimum of n + 1 number [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF] ( θ = min{ θi : 0 ≤ i ≤ n} see Lemma 23 in the Appendix). Furthermore, they proved, that if the matrix M is a P * (κ)-matrix, then θ * p (κ) and θ * c (κ) (defined in the following lemmas) give a feasible predictor and corrector step length pair. The following lemmas show that if either θ * p (κ) or θ * c (κ) is not a feasible step length, than the matrix is not a P * (κ)-matrix. Lemma 15 If there exists an index i (0

≤ i ≤ n) such that θi ≤ θ * p (κ) := 2 (1 -γ)γ (1 + 4κ)n + 2 ,
then matrix M is not a P * (κ)-matrix and the affine Newton direction is a certificate for this. Proof: For any κ ≥ 0 and n ≥ 1

θ * p (κ) < 2 1 + √ 1 + 4κ , therefore if θ0 ≤ θ * p (κ)
, then by the definition of θ0 one has

2 1 + 1 -4e T g/n = θ0 < 2 1 + √ 1 + 4κ
, implying e T g/n < -κ, thus i∈I ∆x i ∆s i < -κnµ = -κ x T s. Therefore, by Lemma 11 the matrix M is not a P * (κ)-matrix and the affine Newton direction ∆x is a certificate for this.

If θi ≤ θ * p (κ), where 0 < i ≤ n, then let consider the following inequality, which was proved by Potra and Liu in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF] on p.158:

(1 -γ)γ + ((1 + 4κ)n + 1) 2 + γ(1 -γ) < (1 + 4κ)n + 2. ( 12
)
Based on the proof of [START_REF] Potra | The Mizuno-Todd-Ye algorithm in a larger neighborhood of the central path[END_REF]Thm 3.3] (see Lemma 24 in the Appendix), Lemma 8 and the definition of t, one has

2 (1 -γ)γ (1 + 4κ)n + 2 = θ * p (κ) ≥ θi ≥ 2 1+ 1+(tγ) -1 (4 g ∞ + 4e T g/n) ≥ 2 1 + 1 + (tγ) -1 (4 g ∞ + 1) = 2 (1-γ)γ (1-γ)γ+ ((1+4κ)n+1)(4 g ∞ +1)+γ(1-γ) . (13)
From inequality ( 13) and ( 12) we get

4 g ∞ + 1 > (1 + 4κ)n + 1. (14)
Since (∆x, ∆s) is a solution of system (3) with a = -xs, and using inequality [START_REF] Kojima | A Unified Approach to Interior Point Algorithms for Linear Complementarity Problems[END_REF] with µ n = x T s, one has

∆x∆s ∞ > (1 + 4κ) 4 x T s = 1 + 4κ 4 a √ xs
I + ∆x i ∆s i ≤ 1 4 μe -xs √ xs 2 = 1 4 μ xs μ - μ xs 2 ≤ 1 4 μ 1 -(1 -t)γ (1 -t)γ n, therefore μ(θ) ≤ 1 + 1 -(1 -t)γ 4(1 -t)γ θ 2 μ. (15)
Since θ * c (κ) is an infeasible step length, there exists an index i such that

x(θ * c (κ)) i s(θ * c (κ)) i < γ μ(θ * c (κ)), namely (1 -θ * c (κ))x i si + θ * c (κ)μ + (θ * c (κ)) 2 ∆x i ∆s i < γ μ(θ * c (κ)).
The predictor point (x, s) ∈ D((1t)γ), so xi si ≥ (1t)γ μ. Furthermore, by inequality (15) one has

(1 -θ * c (κ))(1 -t)γ μ + θ * c (κ)μ + (θ * c (κ)) 2 ∆x i ∆s i < γ 1 + 1 -(1 -t)γ 4(1 -t)γ (θ * c (κ)) 2 μ, which implies (θ * c (κ)) 2 ∆x i ∆s i μ < tγ -θ * c (κ)(1 -(1 -t)γ)+ 1 -(1 -t)γ 4(1 -t) (θ * c (κ)) 2 . ( 16
)
One can check, the following equality by substituting the values of t and θ

* c (κ) 0 ≤ (1 -γ)γ 2 ((1 + 4κ)n + 1) 2 = -tγ + θ * c (κ)(1 -(1 -t)γ)- 1 -(1 -t)γ 4(1 -t)γ (1 + 4κ)n + γ (θ * c (κ)) 2 . Therefore - 1 -(1 -t)γ 4(1 -t)γ (1 + 4κ) n (θ * c (κ)) 2 ≥ tγ -θ * c (κ)(1 -(1 -t)γ) + 1 -(1 -t)γ 4(1 -t) (θ * c (κ)) 2 .
Combining this with inequality ( 16), and then considering the proximity measure estimation in the proof of [START_REF] Potra | The Mizuno-Todd-Ye algorithm in a larger neighborhood of the central path[END_REF]Thm 3.3] (see Lemma 25 in the Appendix), we get

∆x i ∆s i < - 1 -(1 -t)γ 4(1 -t)γ (1 + 4κ) nμ ≤ - (1 + 4κ)μ 4 xs μ - μ xs 2 . (17)
Since (∆x, ∆s) is a solution of system (3) with a = μexs, using inequality [START_REF] Potra | The Mizuno-Todd-Ye algorithm in a larger neighborhood of the central path[END_REF], one get

∆x∆s ∞ > (1 + 4κ)μ 4 xs μ - μ xs 2 = 1 + 4κ 4 a √ xs 2 .
Thus, by the second statement of Lemma 11, the matrix M is not a P * (κ)-matrix and the corrector Newton direction ∆x is a certificate for this.

The following lemma proves, that the predictorcorrector algorithm is well defined. Lemma 17 At each iteration when the value of κ is updated, the new value of θ * p (κ) satisfies the inequality θ ≥ θ * p (κ), and the new point (x(θ * c (κ)), s(θ * c (κ))), determined by the new value of the corrector step size θ * c (κ), is in the D(γ) neighborhood. Proof: In the proof of [START_REF] Potra | The Mizuno-Todd-Ye algorithm in a larger neighborhood of the central path[END_REF]Thm 3.3] (see Lemma 26 in the Appendix) we use the P * (κ) property only for the vector ∆x or ∆x. When parameter κ is updated, then we choose the new value in such a way that the inequality in the definition of P * (κ)-matrices (2) holds for the vectors ∆x and ∆x. Therefore the new value of θ * p (κ) satisfies the inequality θ ≥ θ * p (κ), and the new value of θ * c (κ) determines a point in the D(γ) neighborhood. Now we are ready to state the complexity result for the modified predictor-corrector algorithm for general LCPs in case an initial interior point is available. Theorem 18 Let (x 0 , s 0 ) ∈ F 0 such that (x 0 , s 0 ) ∈ D(γ). Then after at most O (1 + κ)n log (x 0 ) T s 0 ε steps, where κ ≤ κ is the largest value of parameter κ throughout the algorithm, the predictor-corrector algorithm generate a point (x, ŝ), such that xT ŝ ≤ ε and (x, ŝ) ∈ D(γ), or provides a certificate that the matrix is not P * (κ). Proof: We follow the proof of the previous complexity theorem (see Theorem 14). If we take a predictor and a corrector step, then by Theorem 27 and Lemma 17 the decrease of the complementarity gap is at least

3 (1 -γ)γ 2((1 + 4κ)n + 2) x T s n .
This expression is a decreasing function of κ, so at each iteration, when we make a predictor and a corrector step, the complementarity gap decreases at least by

3 (1 -γ)γ 2((1 + 4κ)n + 2)
x T s n .

We take at most as many iterations as in the original predictor-corrector IPM with a P * (κ)-matrix. Thus, referring to the complexity theorem of the original algorithm (see Theorem 28 in the Appendix) we have proved the theorem.

Summary

In this paper we have presented modifications of affine scaling and the predictor-corrector interior point algorithms that enable us to solve general LCPs without the prerequisite to verify special properties of the coefficient matrix. In particular, we have given two constructive proofs of the following EP type theorem from paper [START_REF] Illés | A polynomial path-following interior point algorithms for general linear complementarity problems[END_REF]. We assume that the data are rational (solving problems with computer this is a reasonable assumption), ensuring polynomial encoding size of certificates and polynomial complexity of the algorithms. Theorem 19 Let an arbitrary matrix M ∈ Q n×n , a vector q ∈ Q n and a point (x 0 , s 0 ) ∈ F 0 be given. Then one can verify in polynomial time that at least one of the following statements hold [START_REF] Dikin | Iterative solution of problems of linear and quadratic programming[END_REF] problem LCP has a feasible complementary solution (x, s) whose encoding size is polynomially bounded. ( 2) the matrix M is not in the class of P * (κ) and there is a certificate whose encoding size is polynomially bounded.

2 ,

 2 so by the second statement of Lemma 11 one has M / ∈ P * (κ), and ∆x is a certificate for this. Now let us analyze the corrector step. Lemma 16 If θ * c (κ) := 2γ (1+4κ)n+1 is such a corrector step length that (x(θ * c (κ)), s(θ * c (κ))) / ∈ D(γ), then the matrix M is not a P * (κ)-matrix and the corrector Newton direction is a certificate for this. Proof: Notice that x(θ)s(θ) = (1θ)xs + θ μ e + θ 2 ∆x∆s and μ(θ) = μ + θ 2 ∆x T ∆s n . From Lemma 7 and the proof of [17, Thm 3.3] (see Lemma 25 in the Appendix) we get ∆x T ∆s ≤

Algorithm 2 -Predictor-corrector algorithm

Input:

an upper bound κ > 0 on the value of κ; an accuracy parameter ε > 0; a proximity parameter γ ∈ (0, 1); an initial point (x 0 , s 0 ) ∈ D(γ); begin x := x 0 , s := s 0 , µ := (x 0 ) T s 0 /n, κ := 0; while x T s ≥ ε do Predictor step t =

1-γ (1+4κ)n+1 ; calculate the affine Newton direction (∆x, ∆s) with a = -xs; if (the Newton direction does not exists, or it is not unique) then return the matrix is not P 0 ; [START_REF] Fukuda | EP theorems and linear complementarity problems[END_REF] if (κ(∆x) is not defined) then return the matrix is not P * ;

% see Lemma 10 if (κ(∆x) > κ) then return the matrix is not P * (κ);

Appendix

To make our paper self contained, we include those results from [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF][START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF] that are needed for our developments. All lemmas, theorems are converted to our notations. Lemma 20 (Lemma 4.3 in [12]) Let M be an arbitrary real matrix, δ a (xs) < τ and let (∆x, ∆s) be the affine scaling direction.

Theorem 21 (Theorem 6.1 in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF]) Let M be a P * (κ)matrix, r > 0, τ > 1 and let (∆x, ∆s) be the affine scaling direction. If (x, s) ∈ F 0 , δ a (xs) ≤ τ , 0 ≤ θ and

) is strictly feasible and δ a (x(θ)s(θ)) ≤ τ . Theorem 22 (Corollary 6.1 in [START_REF] Illés | Polynomial affinescaling algorithms for P * (κ) linear complementarity problems[END_REF]) . Lemma 23 (From expressions (3.16), (3.17) in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]) Let M be an arbitrary matrix, (x, s) ∈ D(γ), (∆x, ∆s) be the predictor direction in the predictor-corrector algorithm and let the predictor step length be θ

where

Then we have θ = min θi : 0 ≤ i ≤ n .

Lemma 24 (From the proof of Theorem 3.3 in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]) Let the assumptions of Lemma 23 hold, and θi , 1 ≤ i ≤ n be as it is given in Lemma 23, then

Lemma 25 (From the proof of Theorem 3.3 in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]) Let M be an arbitrary matrix and let the point after the predictor step in the predictor-corrector algorithm satisfy (x, s) ∈ D((1t)γ). Then

Lemma 26 (From Theorem 3.3 in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]) Let M be a P * (κ)-matrix and (x, s) ∈ D(γ). Then the predictor step length satisfy µ g , where µ + g = (x + ) T s + /n = x(θ + ) T s(θ + )/n. Theorem 28 (Corollary 3.4 in [START_REF] Potra | Predictor-corrector methods for sufficient linear complementarity problems in a wide neighborhood of the central path[END_REF]) Let M be a P * (κ)-matrix and (x 0 , s 0 ) be a feasible interior point such that (x 0 , s 0 ) ∈ D(γ). Then in at most O (1 + κ)n log (x 0 ) T s 0 ε steps the predictor-corrector algorithm produces a point (x, ŝ) such that (x, ŝ) ∈ D(γ) and xT ŝ ≤ ε.