Aditya Gahlawat

Giorgio Valmorbida

Stability Analysis of Linear Partial Differential Equations with Generalized Energy Functions

We present a method for the stability analysis of a large class of linear Partial Differential Equations (PDEs) in one spatial dimension. We rely on Lyapunov analysis to establish the exponential stability of the systems under consideration. The proposed test for the verification of the underlying Lyapunov inequalities relies on the existence of solutions of a system of coupled differential equations. We illustrate the application of this method using a PDE actuated by a backstepping computed feedback law. Furthermore, for the case of PDEs defined by polynomial data, we formulate a numerical methodology in the form of a convex optimization problem which can be solved algorithmically. We show the effectiveness of the proposed numerical methodology using examples of different types of PDEs.

I. INTRODUCTION

Various physical quantities pertaining to engineering processes evolve over a spatio-temporal domain. Accurate models for the evolution of these processes are given by Partial Differential Equations (PDEs), a few examples of which may be found in thermonuclear fusion [START_REF] Witrant | A control-oriented model of the current profile in tokamak plasma[END_REF], robotic aircraft [START_REF] Paranjape | PDE boundary control for flexible articulated wings on a robotic aircraft[END_REF], and fluid-solid interactions [START_REF] Bucci | Optimal boundary control with critical penalization for a PDE model of fluid-solid interactions[END_REF]. The study of PDEs from a controls perspective can be broadly classified under two sets of methods, which we refer to as direct and indirect methods. Indirect methods apply standard finite dimensional control theory to PDEs. In such methods, the PDE is approximated by a set of Ordinary Differential Equations (ODEs) [START_REF] El-Farra | Analysis and control of parabolic PDE systems with input constraints[END_REF], [START_REF] Goulart | Global stability analysis of fluid flows using sum-of-squares[END_REF]. These methods are also called early lumping methods. Conversely, we say that a method is direct if it does not approximate the system dynamics, and instead takes into consideration the infinite dimensional nature of the system for the purposes of analysis and control. There has been a plethora of results in the application of direct methods to the analysis and control of PDEs. Backstepping [START_REF] Krstic | Backstepping boundary controllers and observers for the slender Timoshenko beam: Part I -Design[END_REF], [START_REF] Krstic | Adaptive boundary control for unstable parabolic PDEs -Part I: Lyapunov design[END_REF], [START_REF] Smyshlyaev | Adaptive boundary control for unstable parabolic PDEs -Part II: Estimation-based designs[END_REF], [START_REF] Smyshlyaev | Adaptive boundary control for unstable parabolic PDEs -Part III: Output feedback examples with swapping identifiers[END_REF], [START_REF] Meurer | Tracking control for boundary controlled parabolic PDEs with varying parameters: Combining backstepping and differential flatness[END_REF] is one such method for PDE stabilization which relies on the construction of an invertible state transformation that maps the PDE to be controlled to an a priori chosen stable PDE. An example of the direct method based on the Lyapunov approach applied to semilinear parabolic and hyperbolic PDEs can be found in [START_REF] Fridman | An LMI approach to H∞ boundary control of semilinear parabolic and hyperbolic systems[END_REF]. A few more examples of Lyapunov based approaches can be found in [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF], where the authors construct boundary observers for hyperbolic systems as applied to flow control, and [START_REF] Argomedo | A strict control Lyapunov function for a diffusion equation with time-varying distributed coefficients[END_REF], wherein the authors develop distributed controllers for parabolic systems with uncertainties with application to thermonuclear fusion.

A fundamental problem in the study of systems whose dynamics are modeled by PDEs is the stability analysis and the determination of convergence rates to equilibrium solutions. Such tasks can be performed with the Lyapunov's second method [START_REF] Datko | Extending a theorem of A. M. Liapunov to Hilbert[END_REF], [START_REF] Movchan | The direct method of Liapunov in stability problems of elastic systems[END_REF]. In this context, the choice of the class of Lyapunov Functional (LF) candidates is critical and should not introduce conservatism in the analysis. On the other hand, when the parameters of the LFs are unknown variables, it should lead to Lyapunov inequalities that are solvable either Aditya Gahlawat is with the Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, Urbana, IL-61801, USA gahlawat@illinois.edu Giorgio Valmorbida is an associate professor at Laboratoire des Signaux et Systèmes, CentraleSupélec, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 3 Rue Joliot Curie, Gif-sur-Yvette 91192, France giorgio.valmorbida@l2s.centralesupelec.fr numerically or analytically. Recently, a number of direct methods for stability analysis were proposed with their respective numerical formulations [START_REF] Valmorbida | Convex solutions to integral inequalities in two-dimensional domains[END_REF], [START_REF] Papachristodoulou | On the analysis of systems described by classes of partial differential equations[END_REF]. While avoiding the truncation of the PDE dynamics, these methods impose a set of basis functions to obtain a numerical verification of the infinite dimensional inequalities, see for instance [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF], where an orthogonal set of polynomials were used, and [START_REF] Valmorbida | Stability analysis for a class of partial differential equations via semidefinite programming[END_REF], where a standard polynomial basis is used. For timedelay systems, a particular class of infinite-dimensional systems, necessity and sufficiency of a class of LFs, called complete quadratic functionals, has been shown in [START_REF] Gu | Stability of time-delay systems[END_REF]Chapter 5].

A. Contribution

We propose a method to assess the exponential stability of the zero solution of a large class of linear PDEs in one spatial domain. We use a Lyapunov based approach to establish stability of systems under consideration. The chosen LF candidates are composed of integrals with quadratic kernels defined on one and two spatial dimensions. We have previously used such LF candidates for the analysis and control of parabolic PDEs with polynomial data in [START_REF] Gahlawat | A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs[END_REF]. The presented approach extends this previous work by not only considering a much larger class of PDEs (not necessarily defined by polynomial data), but by also formulating an analytical test in the form of a system of coupled differential equations and inequalities for stability analysis. The class of PDEs considered in this work contains parabolic PDEs, hyperbolic PDEs, in-domain and boundary coupled PDEs, PDEs with boundary feedback, and Partial (Integro)-Differential Equations (P(I)DEs). Furthermore, the choice of LFs is an extension of the one in the authors' previous work in [START_REF] Valmorbida | Stability analysis for a class of partial differential equations via semidefinite programming[END_REF], wherein the LFs contained only the integrals defined on one dimensional spatial domains.

The choice of LF candidates leads to the formulation of Lyapunov inequalities in the form of integral inequalities containing one and two dimensional integral terms. Such integral inequalities must be verified on subspaces which are defined by their respective boundary conditions. For this purpose we use Green's theorem and the Fundamental Theorem of Calculus (FTC), which allow us to cast the problem of verification of the Lyapunov inequalities in the form of a system of coupled differential equations. This system contains a non-linear Partial Differential Matrix Equation (PDME), an Ordinary Differential Matrix Inequality (ODMI) and two Linear Matrix Inequalities (LMIs). Thus, the existence of a solution to this system of coupled differential equations implies the stability of the systems considered.

As an illustration of applicability, we show how the proposed methodology can analytically construct an LF certificate for the exponential stability of a backstepping controlled parabolic PDE, further motivating the use of the studied class of LFs. Additional potency of the proposed methodology lies in the ease with which it lends itself to the derivation of a computationally efficient numerical method for stability analysis. To be precise, for the case when the class of PDEs under consideration is defined using polynomial data, we parameterize LFs using sum-of-squares polynomials, which allow us to cast the stability analysis problem as a Semi-Definite Programming (SDP) feasibility problem, which is a convex optimization problem that can be solved efficiently. We study five examples on which we test the effectiveness of the proposed numerical method by computing their parametric stability bounds. A preliminary version of the presented method, addressing only polynomial data for a subclass of systems contained in this paper, has been presented in [START_REF] Gahlawat | A semi-definite programming approach to stability analysis of linear partial differential equations[END_REF].

We consider this work as an initial step towards a unified theory for the analysis and control of linear PDEs akin to the LMI framework for ODEs [START_REF] Dullerud | A course in robust control theory: A convex approach[END_REF]. Indeed, the stability analysis problem examples considered in this paper allow us to validate the choice of the class of LFs, and thus, the proposed framework. Even though the studied numerical examples are restricted to the class of PDEs parametrized by polynomial coefficients, the fact that polynomials can approximate any continuous function with an arbitrary accuracy [START_REF] Kreyszig | Introductory functional analysis with applications[END_REF] in bounded intervals leads us to surmise that the use of polynomials to parameterize the class of LF candidates is valid whenever the actual solutions to the set of conditions are bounded functions.

The presented formulation has interesting connections to results in existing literature. If one were to disregard the double integral terms in the chosen LF candidates, then one recovers the standard LF candidates considered in, for e.g., [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF], [START_REF] Coron | Stabilization of a rotating body beam without damping[END_REF], and [START_REF] Fridman | An LMI approach to H∞ boundary control of semilinear parabolic and hyperbolic systems[END_REF]. The cited results use weighted norms on respective Hilbert spaces as LF candidates to generalize the 'energy' of a system. However, as we show (via numerical experiments) in our previous work in [START_REF] Gahlawat | A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs[END_REF], such LF candidates are inherently conservative. To the best of our knowledge, with the exception of our previous work, the choice of LF candidates, and hence the proposed framework, is entirely new. Even though, to the best of our knowledge, no Lyapunov based framework in the literature considers the class of LF candidates we study, there is an implicit connection between our framework and the backstepping method. For example, in [START_REF] Gahlawat | A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs[END_REF], we show that that the backstepping controlled advection-reaction parabolic PDE admits an LF candidate of the form we consider in this work. This fact is also highlighted by the analytical example considered in Section IV-A. Similarly, it can be shown that various examples of PDEs stabilized by backstepping, as in [START_REF] Bribiesca-Argomedo | Backstepping-forwarding control and observation for hyperbolic PDEs with fredholm integrals[END_REF] and [START_REF] Krstic | Boundary controllers and observers for the linearized Schrödinger equation[END_REF], admit LF candidates of the type we consider.

B. Outline

We begin with Section II in which we state the problem, define the class of PDEs under consideration and derive the integral inequalities which define the Lyapunov inequalities for the stability analysis. Furthermore, we outline the strategy we use in the manuscript to solve the stated problem. In Section III-A we present a method of constructing non-negative and strictly positive integral inequalities on Hilbert Spaces. In Section III-B we use Green's theorem and the FTC to construct slack integrals, that we define to be integral expressions which are identically zero on the subspaces on which the Lyapunov inequalities have to be verified. In Section IV we combine the results from the previous sections to formulate the main contribution of the manuscript. Furthermore, we apply the proposed method to analytically prove the exponentially stability of a parabolic PDE under backstepping control feedback. Finally, in Section V we formulate a numerical methodology for PDEs defined using polynomial data and test the method on various numerical examples.

C. Notation

We denote by

Ω = {(x, y) ∈ R 2 : 0 ≤ y ≤ x ≤ 1}, Ω = {(x, y) ∈ R 2 : 0 ≤ x ≤ y ≤ 1} and Ω = [0, 1] × [0, 1] (Ω = Ω ∪ Ω).
For any γ, n ∈ N and domain Π ⊂ R n , we denote by C γ (Π) the set of γ-times continuously differentiable functions on Π. We also write C 0 (Π) = C(Π). In the following definitions α, β ∈ N. For

w : [0, 1] → R β , w ∈ C α ([0, 1]
), we denote by ∂ i

x w(x) the i-th derivative of w and

wα(x) = w(x) ⊤ ∂xw(x) ⊤ • • • ∂ α x w(x) ⊤ ⊤ , w b α = wα-1(1) ⊤ wα-1(0) ⊤ ⊤ , wα(x) = [wα(x) ⊤ w b α ].
Thus, wα : [0, 1] → R β(α+1) , w b α ∈ R 2βα and wα : [0, 1] → R β(3α+1) . We denote by N ∂ , N ∈ N βα×β(α+1) and N0, N1 ∈ N βα×2βα the matrices satisfying

∂xwα-1(x) = N ∂ wα(x), wα-1(x) = N wα(x), wα-1(0) = N0w b α , wα-1(1) = N1w b α .
We denote

H α [0, 1]; R β = w : [0, 1] → R β : w, ∂xw, . . . , ∂ α-1
x w are absolutely continuous on [0, 1] and

ˆ1 0 (∂ α x w(x)) ⊤ (∂ α x w(x)) dx < ∞ .
We also write

L2 [0, 1]; R β = H 0 [0, 1]; R β .
The space H α [0, 1]; R β is endowed with the inner product and norm

w H α = wα, wα H α = ˆ1 0 wα(x) ⊤ wα(x)dx,
and the space L2 [0, 1]; R β has the norm and inner product

• L 2 = • H 0 and •, • L 2 = •,
• H 0 , respectively. We denote by L ∞ (Φ1, Φ2) the set of Lebesgue measurable and essentially bounded functions mapping Φ1 to Φ2, where Φ1 and Φ2 are any arbitrary sets.

For any m, n ∈ N, we denote by 0m,n the matrix of zeros of dimensions m-by-n and 0m when m = n. Similarly, we denote by In the identity matrix of dimensions n-by-n. For any square matrix Q, we denote He(Q) = 1 2 Q + Q ⊤ . We define the set S n = {M ∈ R n×n : M = M ⊤ } and we say that a matrix valued function S : [0, 1] → S n , is Positive Semi-Definite (PSD) if S(x) 0, for all x ∈ [0, 1].

We denote by S n [x] and S n [(x, y)] the sets of symmetric polynomial matrices of size n-by-n in variables x, and x and y, respectively. Similarly, we denote by R m×n [x] and R m×n [(x, y)] the sets of real polynomial matrices of size m-by-n in variables x and y. We denote by Σ n [x] : [0, 1] → S n the set of Sum-of-Squares (SOS) polynomials in the variable x. Note that, by definition, an SOS polynomial is non-negative for all x ∈ R [5, Chapter 3]. Given α, β, d ∈ N, and the vector of monomials in x and y up to degree d, z ∈ R 1 2 (d+2)(d+1) [(x, y)], (for example, for d = 2, we have z(x, y) = 1 x y x 2 xy y 2 ⊤ ), we define

Z q(α,β,d) (x, y) = I β(α+1) ⊗ z(x, y) ∈ R q(α,β,d)×β(α+1) , (1)
where q(α, β, d) = 1 2 β(α + 1)(d + 2)(d + 1) and ⊗ denotes the Kronecker product.

For any bivariate function K : Ω → R n×n , we define the linear map

Γ [K] = K(x, y), x ≥ y K(y, x) ⊤ , y > x ,
thus satisfying, for any u :

[0, 1] → R n , ˆΩ u(x) ⊤ Γ [K] u(y)dydx = ˆ1 0 ˆx 0 u(x) ⊤ K(x, y)u(y)dydx + ˆ1 0 ˆ1 x u(x) ⊤ K(y, x) ⊤ u(y)dydx.

D. Notational Sets and Maps

To provide a concise presentation, we define the following sets and maps for any α, β, s ∈ N. We begin with the definition of the following sets 1) We say that Z ∈ Π1(s, α, β),

if Z : Ω → R s×β(α+1) and Z ∈ C(Ω). 2) We say that S ∈ Π2(s, α, β),

if S : [0, 1] → S β(3α+1)+2s and is partitioned as

S(x) =     S11(x) S12(x) 0 2βα,s S13(x) 0 2βα,s S12(x) ⊤ 0 s,2βα S22 S23 S13(x) ⊤ 0 s,2βα S ⊤ 23 S33     , for some S11 : [0, 1] → S β(3α+1) , S11 ∈ C([0, 1]), S12, S13 : [0, 1] → R β(α+1)×s , S12, S13 ∈ C([0, 1]), S22, S33 ∈ S s and S23 ∈ R s×s . 3) We say that {K1, K2, H1, H2, B1, B2} ∈ Π3(α, β) (4)
if

K1 : [0, 1] → S βα , K2 : [0, 1] → R βα×2βα , K1, K2 ∈ C 1 ([0, 1]), H1, H2 : Ω → R βα×βα , H1, H2 ∈ C 1 (Ω), B1 : [0, 1] → R β(α+1)×βα , B1 ∈ C([0, 1]), B2 ∈ R 2βα×βα .
4) For each of the sets Πi, i ∈ {1, 2, 3}, we write Π R i when the functions under consideration are polynomial in their respective arguments. With the set definitions provided, we now define the following maps:

1) Given any α, β ∈ N, positive scalar δ, and

J b : [0, 1] → S β , J : Ω → R β×β , L1 : [0, 1] → R β×β(α+1) , L2 : Ω → R β×β(α+1) , L3 : Ω → R β×β(α+1) ,
we say that

{M b , M } = Ξ1(L1, L2, L3, J b , J , α, β, δ) (5)
if

M b (x) = -He J b (x)L1(x) 0 β,2βα 0 3βα,β(α+1) 0 3βα,2βα -He δJ b (x) 0 β,3βα 0 3βα,β 0 3βα , M (x, y) = - 1 2 J(x, y) ⊤ L1(x) 0 βα,β(α+1 ) ⊤ + J (x, y)L1(y) 0 βα,β(α+1) - 1 2 J b (x)L2(x, y) 0 βα,β(α+1) + J b (y)L3(y, x) 0 βα,β(α+1) ⊤ - 1 2 M (x, y) + 2δ J (x, y) 0 β,βα 0 βα,β 0 βα , with M (x, y) = ˆy 0 J(x, z) 0 βα,β L3(z, y) + J (y, z) 0 βα,β L3(z, x) ⊤ dz + ˆx y J (x, z) 0 βα,β L2(z, y) + J (z, y) ⊤ 0 βα,β L3(z, x) ⊤ dz + ˆ1 x J(z, x) ⊤ 0 βα,β L2(z, y) + J (z, y) ⊤ 0 βα,β L2(z, x) ⊤ dz.
2) Given any Z ∈ Π1(s, α, β) and S ∈ Π2(s, α, β) we denote by

{S b , S} = Ξ2(Z, S) (6)
if

S b (x) =S11(x), S(x, y) =S12(x)Z(x, y) + Z(y, x) ⊤ S13(y) ⊤ + ˆy 0 Z(z, x) ⊤ S33Z(z, y)dz + ˆx y Z(z, x) ⊤ S ⊤ 23 Z(z, y)dz + ˆ1 x Z(z, x) ⊤ S22Z(z, y)dz.
3) Given any α, β ∈ N, {K1, K2, H1, H2, B1, B2} ∈ Π3(α, β) and

O1 : [0, 1] → R βα×β(α+1) , O2 ∈ R βα×2βα , we say that {K b , H b , H, B b , B} = Ξ3(K1, K2, H1, H2, B1, B2, O1, O2), (7) if
K b = K b1 (x) K b2 (x) 0 2βα,β(α+1) K b3 , H b (x) = H b1 (x) H b2 (x) H b3 (x) 0 2βα , H(x, y) = 1 2 N ⊤ (∂yH1(x, y) -∂xH2(x, y)) N + 1 2 N ⊤ H1(x, y)N ∂ -N ⊤ ∂ H2(x, y)N , B b (x) = 0 β(α+1) B1(x)O2 B2O1(x) B2O2 , B(x, y) = 1 2 B1(x)O1(y) + O1(x) ⊤ B1(y) ⊤ ,
where

K b1 (x) = N ⊤ ∂xK1(x)N + N ⊤ K1(x)N ∂ + N ⊤ ∂ K1(x)N, K b2 (x) = N ⊤ ∂ K2(x) + N ⊤ ∂xK2(x), K b3 = N ⊤ 0 K1(0)N0 -N ⊤ 1 K1(1)N1 + N0K2(0) -N1K2(1), H b1 (x) = -N ⊤ (H1(x, x) + H2(x, x)) N, H b2 (x) = N ⊤ H1(x, 0)N0, H b3 (x) = N ⊤ 1 H2(1, x)N,
and where the matrices N , N0, N1, and N ∂ are defined in Section I-C.

II. PROBLEM STATEMENT

We study the following class of linear Partial Differential Equations (PDEs) ∂tw(x, t) = A1(x)wα(x, t) + ˆx 0 A2(x, y)wα(y, t)dy

+ ˆ1 x A3(x, y)wα(y, t)dy, (8a)
w(•, t) ∈ B, ∀t ≥ 0, (8b)
where α+1) . The boundary conditions define the set

A1 : [0, 1] → R β×β(α+1) , A2 : Ω → R β×β(α+1) , A3 : Ω → R β×β(
B = u ∈ H α [0, 1]; R β : ˆ1 0 F1(x) F2 ūα(x)dx = 0 βα,1 , (8c)
where F1 : [0, 1] → R βα×β(α+1) and F2 ∈ R βα×2βα . If the terms A2 and A3 are identically zero, then one recovers the standard class of PDEs with local dynamics as, for instance, the reactionadvection-diffusion and beam equation [START_REF] Evans | Partial differential equations[END_REF]. However, the integral terms allow us to consider PDEs with non-local behavior. In fact, such PDEs model various processes in engineering and biology like microelectro-mechanical-systems, Ohmic heating phenomena, chemotaxis, and cell dynamics [START_REF] Kavallaris | Non-Local Partial Differential Equations for Engineering and Biology[END_REF].

The motivation to study PDEs of the form in ( 8) lies in the fact that various types of PDEs belong to this class. Each of the following PDEs can be cast in the form of [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF] and is parameterized by a positive scalar λ which can alter the respective stability property. The parameter λ will later help us analyze the effectiveness of the methodology we develop in the paper.

Example 1: Partial (Integro) Differential Equation (P(I)DE): We begin with the following P(I)DE

∂tv(x, t) =∂xv(x, t) + λ ˆx 0 (x -y)v(y, t)dy + λ ˆ1 x (x + y)v(y, t)dy, (9a)
v(1, t) =0. (9b)
Using a strongly continuous semigroup approach as in [32, Theorem 2.1] and the perturbation result in [START_REF] Curtain | An introduction to infinite-dimensional linear systems theory[END_REF]Theorem 3.2.1], it can be shown that (9) admits a unique classical solution for a sufficiently regular initial condition. This example simplifies the system studied in [START_REF] Bribiesca-Argomedo | Backstepping-forwarding control and observation for hyperbolic PDEs with fredholm integrals[END_REF], where a strategy for boundary control is presented.

Example 2: Wave Equation: Consider the following hyperbolic equation

∂ 2 t v(x, t) = ∂ 2 x v(x, t), (10a) ∂xv(0, t) -(1 -λ)∂tv(0, t) = 0, (10b) v(1, t) = 0. (10c)
Using separation of variables, it can be shown that the analytical solution to the PDE is given by the following convergent series

v(x, t) = ∞ n=1 e σnt φn(x)an, (11)
where 2+x) , and an are scalars dependent on the initial conditions.

Re(σn) = -1 2 ln 2-λ λ , φn(x) = e σnx -e -σn(
Example 3: Parabolic PDE with Scalar Coefficients: We now consider the following parabolic PDE

∂tv(x, t) = ∂ 2 x v(x, t) + λv(x, t), v(0, t) = 0, ∂xv(1, t) = 0. ( 12
)
Using separation of variables we calculate the solution to this PDE as

v(x, t) = ∞ n=1 e σnt φn(x)an, (13)
where

σn = λ -(2n -1)π 2 /4, φn = √ 2 sin ((2n -1)πx/2
) and the scalar coefficients an depend on the initial condition.

(x, t) = -j∂ 2 x v(x, t), (14a) ∂xv(0, t) -j(1 -λ)v(0, t) = 0, v(1, t) = 0, (14b)
where

v : [0, 1] × [0, ∞) → C
(x, t) = ∂x [φ(x)∂xv(x, t)] + θ(x)v(x, t), (15a)
∂xv(0, t) + ˆ1 0 ψ(x)v(x, t)dx = 0, (15b) v(1, t) + ∂xv(1, t) + ˆ1 0 ξ(x)v(x, t)dx = 0, (15c)
where

φ(x) = x 2 + 5x + 1, θ(x) = λ -x, ψ(x) = x 2 + 1 and ξ(x) = x.
The well-posedness of this equation can be established using arguments as in [START_REF] Balogh | Stability of partial difference equations governing control gains in infinite-dimensional backstepping[END_REF]Section 6]. This example allows to illustrate the applicability of the numerical methods developed in this paper to systems with polynomial spatially-distributed coefficients and boundary feedback.

We would like to remark that deriving general conditions on A1, A2, A3, F1, and F2 which render (8) well-posed is beyond the scope of the paper. Instead, as for each of the examples above, the well-posedness needs to be established individually for each type of equation considered. Thus, for the remainder of the exposition, we assume that ( 8) is well posed, that is, the PDE admits a unique classical solution which is continuously differentiable in time and belongs to H α [0, 1]; R β for all t > 0.

In this work we wish to establish the stability of (8) by constructing Lyapunov Functional (LF) certificates of exponential stability. In particular, we consider the following class of functions as LF candidates

V (u) := 1 2 ˆ1 0 u(x) ⊤ T b (x)u(x)dx + 1 2 ˆΩ u(x) ⊤ Γ T u(y)dydx, u ∈ L2 [0, 1], R β , (16)
where T b : [0, 1] → S β and T : Ω → R β×β are Lebesgue measurable and satisfy

T b ∈ L ∞ [0, 1], S β , T ∈ L ∞ Ω, R β×β .
In order to formulate the conditions for the exponential stability of (8), we define the following integral expression. For any scalar δ > 0, we define

V d (u) = ˆ1 0 ūα(x) ⊤ U b (x)ūα(x)dx + ˆΩ uα(x) ⊤ Γ Ū uα(y)dydx, ∀u ∈ H α [0, 1], R β , (17)
where

{U b , Ū } = Ξ1(A1, A2, A3, T b , T , α, β, δ),
and where the map Ξ1 is defined in [START_REF] Blekherman | Semidefinite optimization and convex algebraic geometry[END_REF].

In the Appendix B we show that for solutions w of the PDE (8) we have

V d (w(•, t)) = -∂tV (w(•, t)) -2δV (w(•, t)).
We now state the conditions for the exponential stability of (8) in the following lemma.

Lemma 1. Given the PDE [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF], suppose there exist positive scalars µ1, µ2, δ and Lebesgue measurable functions

T b : [0, 1] → S β and T : Ω → R β×β , T b ∈ L ∞ [0, 1], S β , T ∈ L ∞ Ω, R β×β such that µ1 u 2 L 2 ≤ V (u) ≤ µ2 u 2 L 2 , ∀u ∈ L2 [0, 1]; R β , (18a) V d (u) ≥ 0, ∀u ∈ B ⊂ H α [0, 1]; R β , ( 18b
)
where V and V d are defined in [START_REF] Evans | Partial differential equations[END_REF] and [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF], respectively. Then, there exists a positive scalar κ such that the solutions of (8) satisfy

w(•, t) L 2 ≤ κe -δt w(•, 0) L 2 , ∀t ≥ 0. ( 19
)
Note that the above result gives the exponential decay of the L2norm of the state. Without using the bounds in terms of the L2 norm as in (18a), one can still obtain the decay of the form V(w(•, t)) ≤ e -2δt V(w(•, 0)), for all t ≥ 0.

We have reduced the problem of assessing the exponential stability to the problem of verification of the integral inequalities in [START_REF] Fridman | An LMI approach to H∞ boundary control of semilinear parabolic and hyperbolic systems[END_REF]. Verification of (18a) is straightforward. We generalize the results in [START_REF] Gahlawat | A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs[END_REF] and [START_REF] Valmorbida | Stability analysis for a class of partial differential equations via semidefinite programming[END_REF] to verify

µ1 u 2 L 2 ≤ V (u) ≤ µ2 u 2 L 2 , ∀u ∈ L2 [0, 1]; R β ,
using a PSD matrix valued function, which ensures that the integral inequalities hold. Verification of (18b) is unfortunately not as straight forward. One method is to follow a similar procedure as the one adopted for (18a). That is, if we can construct an integral expression R (u) which satisfies

R (u) ≥ 0, ∀u ∈ H α [0, 1]; R β , (20)
to be a lower bound satisfying

V d (u) -R (u) ≥ 0, ∀u ∈ H α [0, 1]; R β ,
then, the existence of such functional R (u) would therefore imply

V d (u) ≥ 0, ∀u ∈ H α [0, 1]; R β , (21)
thus ensuring that (18b) holds since B is a subset of H α [0, 1]; R β . However, this approach would be conservative since (18b) requires V d (u) ≥ 0 to hold on the subset B ⊂ H α [0, 1]; R β , and not on H α [0, 1]; R β as in [START_REF] Goulart | Global stability analysis of fluid flows using sum-of-squares[END_REF]. To solve this problem, we rely on Slack Integrals, S(u), which are integral expressions that satisfy

S(u) = 0, ∀u ∈ B. (22)
Then, to verify (18b) we may construct the term R (u) as in [START_REF] Gahlawat | A semi-definite programming approach to stability analysis of linear partial differential equations[END_REF] and a slack integral S(u) such that

V d (u) + S(u) -R (u) ≥ 0, ∀u ∈ H α [0, 1]; R β . ( 23
)
If this condition holds, then, owing to [START_REF] Gahlawat | A semi-definite programming approach to stability analysis of linear partial differential equations[END_REF],

V d (u) + S(u) ≥ 0, ∀u ∈ H α [0, 1]; R β ,
and since the slack integral S(u) satisfies [START_REF] Gu | Stability of time-delay systems[END_REF], we obtain

V d (u) ≥ 0, ∀u ∈ B,
thus verifying (18b). To summarize, it is the slack integrals S(u) that allow us to additionally consider the effects of the boundary conditions u(t) ∈ B, thus making it possible to verify that V d (u) ≥ 0 for systems for which damping is introduced by the boundaries. We use quadratic forms of the FTC and Green's theorem to construct slack integrals. Such a formulation allows us to reduce the verification of (23) to a system of coupled differential equations containing a Partial Differential Matrix Equation (PDME) on Ω, an Ordinary Differential Matrix Inequality (ODMI) on the interval [0, 1] and two algebraic Linear Matrix Inequalities (LMIs). Thus, the existence of a solution to the system of coupled differential equations would ensure that (18b) holds. The construction of positive/nonnegative integral inequalities is presented in Section III-A. The formulation of slack integrals is presented in Section III-B. Finally, in Section IV we formulate the conditions to verify [START_REF] Fridman | An LMI approach to H∞ boundary control of semilinear parabolic and hyperbolic systems[END_REF].

III. PRELIMINARIES

In this section we construct integral inequalities on Hilbert spaces and slack integrals which will help us in verifying [START_REF] Fridman | An LMI approach to H∞ boundary control of semilinear parabolic and hyperbolic systems[END_REF].

A. Integral Inequalities on Hilbert Spaces

We begin by constructing non-negative integral terms on

H α [0, 1]; R β . Lemma 2. Given any α, β, p ∈ N, Yp ∈ Π1(p, α, β), and R ∈ Π2(p, α, β), consider R (u) = ˆ1 0 ūα(x) ⊤ R b (x)ūα(x)dx + ˆΩ uα(x) ⊤ Γ R uα(y)dydx, u ∈ H α [0, 1]; R β , ( 24
)
where {R b , R} = Ξ2(Yp, R) and Π1, Π2, Ξ2 are defined in (2), ( 3), [START_REF] Bribiesca-Argomedo | Backstepping-forwarding control and observation for hyperbolic PDEs with fredholm integrals[END_REF], respectively. If

R(x) 0, ∀x ∈ [0, 1], ( 25
) then R (u) ≥ 0, ∀u ∈ H α [0, 1]; R β . (26)
The proof of this lemma is provided in Appendix B. We now present a similar result in which we construct integral inequalities on L2 [0, 1]; R β . Proposition 1. Given any β, r ∈ N, Vr ∈ Π1(r, 0, β), and T ∈ Π2(r, 0, β), consider

T (u) = ˆ1 0 u(x) ⊤ T b (x)u(x)dx + ˆΩ u(x) ⊤ Γ T u(y)dydx, u ∈ L2 [0, 1]; R β , (27)
where {T b , T } = Ξ2(Vr, T ). Suppose that T (x) satisfies one of the following conditions: 1) there exists a positive scalar ǫ such that

T (x) - ǫI β 0 β,2r 0 2r,β 02r 0, ∀x ∈ [0, 1]; (28a)
2) there exists a positive scalar ǫ, matrix valued function Q1 :

[0, 1] → S β , Q1 ∈ C([0, 1
]), and matrix Q2 ∈ R β×r such that

T (x) =   Q1(x) Q ⊤ 2 0 r,β     Q1(x) Q ⊤ 2 0 r,β   ⊤ , Q1(x) ǫI β , (28b)
for all x ∈ [0, 1]; 3) there exists a positive scalar ǫ, matrix valued function Q1 :

[0, 1] → S β , Q1 ∈ C([0, 1]), and matrix Q2 ∈ R β×r such that

T (x) =   Q1(x) 0 r,β Q ⊤ 2     Q1(x) 0 r,β Q ⊤ 2   ⊤ , Q1(x) ǫI β , (28c)
for all x ∈ [0, 1]; then, there exist positive scalars θ1, θ2 such that

θ1 u 2 L 2 ≤ T (u) ≤ θ2 u 2 L 2 , ∀u ∈ L2 [0, 1]; R β .
The proof of this proposition can be found in Appendix B.

B. Slack Integrals

In this section we construct slack integrals, which we define as integral expressions S(u) satisfying S(u) = 0, ∀u ∈ B, with B ⊂ H α [0, 1]; R β given in (8c). We formulate the slack integrals using the FTC, Green's theorem, and the definition of the boundary conditions in (8c).

The following results use the set Π3 and the map Ξ3 defined in ( 4) and [START_REF] Bucci | Optimal boundary control with critical penalization for a PDE model of fluid-solid interactions[END_REF], respectively. Lemma 3 (FTC quadratic form). Given any α, β ∈ N and

{K 1 , K 2 , •, •, •, •} ∈ Π 3 (α, β), the following identity holds ˆ1 0 ūα(x) ⊤ He (K b (x)) ūα(x)dx = 0, ∀u ∈ H α [0, 1]; R β , ( 29
)
where {K b , •, •, •, •} = Ξ 3 (K 1 , K 2 , •, •, •, •, •, •).
Proof: The identity ( 29) is established by developing

ˆ1 0 d dx l(x) dx -(l(1) -l(0)) = 0, with l(x) =uα-1(x) ⊤ K1(x)uα-1(x) + uα-1(x) ⊤ K2(x)u b α .
Next, we present the quadratic form of the Green's theorem. The proof of the following lemma is provided in the Appendix B.

Lemma 4 (Green's theorem quadratic form). Given any α, β ∈ N and {•, •, H1, H2, •, •} ∈ Π3(α, β), the following identity holds

ˆ1 0 ūα(x) ⊤ He (H b (x)) ūα(x)dx + ˆΩ uα(x) ⊤ Γ H uα(y)dydx = 0, ∀u ∈ H α [0, 1]; R β , where {•, H b , H, •, •} = Ξ3(•, •, H1, H2, •, •, •, •).
In the following lemma we formulate an integral equation that holds on the set B defined in (8c), the proof of which is provided in Appendix B. Lemma 5. Consider the set B defined by F1 : [0, 1] → R βα×β(α+1) and F2 ∈ R βα×2βα as in (8c). For any {•, •, •, •, B1, B2} ∈ Π3(α, β) the following identity holds true for all u ∈ B

ˆ1 0 ūα(x) ⊤ He (B b (x)) ūα(x)dx + ˆΩ uα(x) ⊤ Γ B uα(y)dydx = 0, where {•, •, •, B b , B} = Ξ3(•, •, •, •, B1, B2, F1, F2).
We now present the main result of this section wherein we use the results in Lemmas 3-5 to construct slack integrals.

Theorem 1. Consider the set B defined by F1 : [0, 1] → R βα×β(α+1) and F2 ∈ R βα×2βα as in (8c). For any {K1, K2, H1, H2, B1, B2} ∈ Π3(α, β), define

S(u) = ˆ1 0 ūα(x) ⊤ He (K b (x) + H b (x) + B b (x)) ūα(x)dx + ˆΩ uα(x) ⊤ Γ H + B uα(y)dydx, ( 30
)
where {K b , H b , H, B b , B} = Ξ3(K1, K2, H1, H2, B1, B2, F1, F2).

The following identity holds true

S(u) = 0, ∀u ∈ B. ( 31
)
Proof: We begin by considering the following decomposition

S(u) = ˆ1 0 ūα(x) ⊤ He (K b (x) + H b (x) + B b (x)) ūα(x)dx + ˆΩ uα(x) ⊤ Γ H + B uα(y)dydx = 3 i=1 Θi, (32)
where

Θ1 = ˆ1 0 ūα(x) ⊤ He (K b (x)) ūα(x)dx, Θ2 = ˆ1 0 ūα(x) ⊤ He (H b (x)) ūα(x)dx + ˆΩ uα(x) ⊤ Γ H uα(y)dydx, Θ3 = ˆ1 0 ūα(x) ⊤ He (B b (x)) ūα(x)dx + ˆΩ uα(x) ⊤ Γ B uα(y)dydx.
From Lemmas 3 and 4 we have that

Θ1 = 0 and Θ2 = 0, ∀u ∈ H α [0, 1]; R β . ( 33
)
From Lemma 5 we have that

Θ3 = 0, ∀u ∈ B. (34)
Therefore, from ( 32)- [START_REF] Paranjape | PDE boundary control for flexible articulated wings on a robotic aircraft[END_REF] we conclude that the expression in [START_REF] Movchan | The direct method of Liapunov in stability problems of elastic systems[END_REF] holds for all u ∈ B.

IV. MAIN RESULT

In this section we use the results of the previous two sections to formulate the conditions for exponential stability of [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF]. We present the following theorem.

Theorem 2. Consider the PDE [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF]. Given any positive scalars ǫ, δ and p, r ∈ N, suppose there exist functions and matrices Vr ∈ Π1(r, 0, β), T ∈ Π2(r, 0, β), (35a)

Yp ∈ Π1(p, α, β), R ∈ Π2(p, α, β), (35b) {K1, K2, H1, H2, B1, B2} ∈ Π3(α, β), (35c)
such that

T (x) satisfies either (28a), (28b), or (28c), (36a) R(x) 0, ∀x ∈ [0, 1], (36b) He (U b (x) + K b (x) + H b (x) + B b (x)) -R b (x) 0, ∀x ∈ [0, 1], (36c) Ū (x, y) + H(x, y) + B(x, y) -R(x, y) = 0α+1, ∀(x, y) ∈ Ω, ( 36d
)
where

{U b , Ū } =Ξ1(A1, A2, A3, T b , T , δ), {T b , T } =Ξ2(Vr, T ), {R b , R} =Ξ2(Yp, R), {K b , H b , H, B b , B} =Ξ3(K1, K2, H1, H2, B1, B2),
and where Πi, Ξi, i ∈ {1, 2, 3}, are defined in (2)- [START_REF] Bucci | Optimal boundary control with critical penalization for a PDE model of fluid-solid interactions[END_REF].

Then, PDE (8) is exponentially stable.

Proof: Since (36a) holds, we conclude from Proposition 1 that there exist positive scalars θ1, θ2 such that

θ1 u 2 L 2 ≤ T (u) ≤ θ2 u 2 L 2 , ∀u ∈ L2 [0, 1]; R β , (37)
where T (u) is defined by T b and T as in [START_REF] Krstic | Adaptive boundary control for unstable parabolic PDEs -Part I: Lyapunov design[END_REF]. Now, let us define

P (u) = ˆ1 0 ūα(x) ⊤ P b (x)ūα(x)dx + ˆΩ uα(x) ⊤ Γ P uα(y)dydx, with P b (x) =He (U b (x) + K b (x) + H b (x) + B b (x)) -R b (x), P (x, y) = Ū (x, y) + H(x, y) + B(x, y) -R(x, y).
Then, as a consequence of (36c)-(36d) we have

P (u) ≥ 0, ∀u ∈ H α [0, 1]; R β . (38)
Additionally, by definition of P (u) we have

P (u) = V d (u) + S(u) -R (u), (39)
with V d (u) defined as in [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF], S(u) defined as in [START_REF] Meurer | Tracking control for boundary controlled parabolic PDEs with varying parameters: Combining backstepping and differential flatness[END_REF], and R (u) defined as in [START_REF] Kreyszig | Introductory functional analysis with applications[END_REF]. Therefore, we may use ( 38) and ( 39) to conclude

V d (u) + S(u) ≥ R (u), ∀u ∈ H α [0, 1]; R β . (40)
Since (36b) holds, from Lemma 2 we conclude

R (u) ≥ 0, ∀u ∈ H α [0, 1]; R β .
Thus, (40) can be reduced to

V d (u) + S(u) ≥ 0, ∀u ∈ H α [0, 1]; R β . (41)
From Theorem 1 we have

S(u) = 0, ∀u ∈ B, therefore (41)
implies

V d (u) ≥ 0, ∀u ∈ B. (42)
Since ( 37) and (42) hold, using Lemma 1 proves that PDE ( 8) is exponentially stable. The inequalities in (36a)-(36b) are Linear Matrix Inequalities (LMIs) for matrix valued functions. The equation in (36d) is a nonlinear Partial Differential Matrix Equation (PDME) since it contains the partial derivatives of the variables H1 and H2 owing to Lemma 4, the non-linear terms in the variables appearing in R, namely products between Yp and blocks of R(x), and the term T , namely products between Vr and blocks of T (x), affecting Ū . Moreover, (36c) is an Ordinary Differential Matrix Inequality (ODMI) since it contains derivatives of variables K1 and K2 due to Lemma 3 and also contains the values of the variables H1 and H2 on the boundary of Ω ⊂ R 2 . Therefore, the conditions of Theorem 2 require the solution of a nonlinear PDME subjected to an ODMI and two LMIs. It is beyond the scope of this paper to establish the well-posedness of the PDME. However, we show in the following subsection that a solution exists for a particular example of a PDE. Moreover, in Section V we show that, for a particular subclass of (8), the verification of the conditions of Theorem 2 can be reduced to a convex feasibility problem, and thus, can be solved numerically.

A. Boundary Controlled Parabolic PDE-Backstepping feedback law

We now show that the conditions of Theorem 2 are verified for a boundary controlled PDE. Consider the following parabolic PDE with a backstepping control feedback law ∂tw(x, t) =∂ 2

x w(x, t) + λw(x, t), (43a)

w(0, t) =0, w(1, t) - ˆ1 0 J(1, x)w(x, t)dx = 0, (43b)
where λ ∈ R, the function J(1, x) is obtained from the backstepping control kernel,

J(x, y) = -(λ + κ)y I1 (λ + κ) (x 2 -y 2 ) (λ + κ) (x 2 -y 2 ) , (44)
where κ ∈ R is any scalar satisfying 0 < κ < ∞ and I1 is the first order modified Bessel function of the first kind. The function J(x, y), as explained in [28, Section 4.7], is the solution to the following kernel PDE

∂ 2 x J(x, y) -∂ 2 y J(x, y) = (λ + κ) J(x, y), (45a)
J(x, 0) = 0, J(x, x) = - 1 2 x (λ + κ) . (45b)
The system ( 43) is in the form of system ( 8) with α = 2, β = 1 and

A1(x) = λ 0 1 , A2(x, y) = A3(x, y) = 0 β,β(α+1) , (46a)
F1(x) = -J(1, x) 0 0 0 0 0 , F2 = 1 0 0 0 0 0 1 0 . ( 46b
)
Remark. The PDE (43) is exponentially stable. Indeed, as explained in [START_REF] Krstic | Boundary control of PDEs: A course on backstepping designs[END_REF]Section 4.7], with the following invertible variable transformation

v(x, t) = (Gw(•, t)) (x) = w(x, t) - ˆx 0 J(x, y)w(y, t)dy, (47)
we obtain

∂tv(x, t) =∂ 2 x v(x, t) -κv(x, t), (48a) v(0, t) =0, v(1, t) = 0, (48b)
which, since κ > 0, is exponentially stable. Therefore, for the particular example considered in this section, the exponential stability in the L2-norm can be proven using a simple energy function applied to the system in the target coordinates (48).

We would like to use Theorem 2 to prove the exponential stability of (43) in the original coordinates (48). That is, by showing the existence of functions [START_REF] Renardy | An introduction to partial differential equations[END_REF] that satisfy [START_REF] Smyshlyaev | Adaptive boundary control for unstable parabolic PDEs -Part II: Estimation-based designs[END_REF].

Set

ǫ ∈ (0, 1), δ = κ, (49a)
T (x) =   1 -1 0 -1 1 0 0 0 0   , Vr(x, y) = J(x, y), r = 1, (49b) R(x) =     R11(x) R12(x) 04,1 R13(x) 04,1 R12(x) ⊤ 01,4 R22 R23 R13(x) ⊤ 01,4 R ⊤ 23 R33     , (49c) R11(x) =     J(x, x) 2 -J(x, x) 0 -J(x, x) 1 0 0 0 0 03,4 04,3 04
    , (49d)
R12(x) =   -J(x, x) 1 0   , R13(x) = 03,1, (49e)
R22 = 1, R23 = R33 = 0, (49f) Yp(x, y) = -∂xJ(x, y) 0 0 , p = 1, (49g)
K1(x) = -J(x, x) 1 2 1 2 0 , K2(x) = 02,4, (49h)
H1(x, y) = ∂yh1(x, y) -h1(x, y) 0 0 , (49i) H2(x, y) = ∂xh2(x, y) 0 h1(x, y) 0 , (49j)
B1(x) = -   [∂xJ(x, y)] x=1 y=x [∂yh1(x, y)] y=0 0 0 0 0   , (49k) B2 =     -J(1, 1) 0 1 0 0 0 0 -1     , (49l)
where h1(x, y) =J(x, y) -ˆ1 x J(z, x)J(z, y)dz, (50)

h2(x, y) =J(x, y) + ˆ1 x J(z, x)J(z, y)dz. ( 51
)
With the above values, we state the following proposition, which, using Theorem 2 allows us to conclude the exponential stability of ( 43). [START_REF] Smyshlyaev | Adaptive boundary control for unstable parabolic PDEs -Part II: Estimation-based designs[END_REF] holds with A1(x), A2(x, y), A3(x, y), F1(x) and F2 as in (46) and the functions in (49).

Proposition 2. The set of equations and inequalities in

Proof: If we define Q1(x) = 1 and Q2 = -1, then it is easily established that T ∈ Π2(r, 0, β) and satisfies (28b) for any ǫ ∈ (0, 1). Thus, (36a) is satisfied.

Similarly, R ∈ Π2(p, α, β) and satisfies [START_REF] Krstic | Boundary controllers and observers for the linearized Schrödinger equation[END_REF] since

R(x) = r(x) ⊤ r(x), r(x) = -J(x, x) 1 01,5 1 0 .
Thus, (36b) is satisfied. Using (49b) and the definition {T b , T } = Ξ1(Vr, T ), we calculate

T (x, y) = -h1(x, y), ( 52
)
where h1 is defined in (50). Using these functions and the definitions in (46), we may apply [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF] to compute

He (U b (x)) =     -(λ + δ) 0 -1 2 0 0 0 -1 2 0 0 03,4 04,3 04
    , (53)
Ū (x, y) =   (λ + δ)h1(x, y) 0 1 2 h1(x, y) 0 0 0 1 2 h1(x, y) 0 0   . ( 54
)
Using the definitions in (49h) and Lemma 3 we calculate

He (K b (x)) =   K b1 (x) 03,2 03,2 02,3 K b2 02 02,3 02 K b3   , (55)
where

K b1 (x) =   -d dx J(x, x) -J(x, x) 1 2 -J(x, x) 1 0 1 2 0 0   , K b2 = J(1, 1) -1 2 -1 2 0 , K b3 = -J(0, 0) 1 2 1 2 0 .
Using the definitions in (49i)-(49j) and Lemma 4 we obtain

He (H b (x)) =     J(x, x) 2 -d dx J(x, x) 0 0 m(x) 0 0 m(x) ⊤ 06     , (56) H(x, y) =   1 2 ∂ 2 y h1 -∂ 2 x h2 0 -1 2 h1(x, y) -∂xJ(x, y) 0 0 -1 2 h1(x, y) 0 0   , (57)
where h1 and h2 are defined in ( 50) and ( 51), respectively, and

m(x) ⊤ =      1 2 [∂xJ(x, y)] x=1 y=x -J(1, 1)J(1, x) 1 2 J(1, x) 1 2 [∂yh1(x, y)] y=0 0     
.

Similarly, we use Lemma 5 for the definitions in (49k)-(49l) to obtain

He (B b (x)) = 03 B b1 (x) B b1 (x) ⊤ B b2 , (58)
B(x, y) =   b(x, y) 0 0 0 0 0 0 0 0   , (59)
where

B b1 (x) = m(x) 02,4 , B b2 =     -J(1, 1) 1 2 0 0 1 2 0 0 0 0 0 0 -1 2 0 0 -1 2 0     , b(x, y) = 1 2 J(1, x) [∂xJ(x, y)] x=1 + 1 2 J(1, y) [∂xJ(x, y)] x=1 y=x .
Now, using ( 53), ( 55), ( 56) and ( 58) we get

He (U b (x) + K b (x) + H b (x) + B b (x)) =     J(x, x) 2 + κ -δ -J(x, x) 0 -J(x, x) 1 0 0 0 0 03,4 04,3 04
    , (60)
where we have used the property in (45b) that

d dx J(x, x) = - 1 2 (λ + κ).
Similarly, using (54), ( 57) and (59) we get

Ū (x, y) + H(x, y) + B(x, y) =   (λ + δ)h1 + 1 2 ∂ 2 y h1 -∂ 2 x h2 + b(x, y) 0 0 -∂xJ(x, y) 0 0 0 0 0   . (61)
Now, let us define

g(x, y) = ˆ1 x ∂zJ(z, x)∂zJ(z, y)dz + J(x, x)∂xJ(x, y). ( 62
)
Then, applying Lemma A.3 in Appendix A produces

Ū (x, y) + H(x, y) + B(x, y) =   (δ -κ)h1(x, y) + g(x, y) 0 0 -∂xJ(x, y) 0 0 0 0 0   . ( 63
)
Using the definitions in (49a), (49c)-(49g) and ( 24) we calculate

R(x, y) =   g(x, y) 0 0 -∂xJ(x, y) 0 0 0 0 0   . ( 64
)
Now, from (49a), (49d) and (60) we get

He (U b (x) + K b (x) + H b (x) + B b (x)) -R b (x) = 07.
Thus, (36c) is satisfied. Similarly, using ( 63) and (64) we get

Ū (x, y) + H(x, y) + B(x, y) -R(x, y) = 03.
Therefore, (36d) holds.

In conclusion, we have proved that for the PDE [START_REF] Yamashita | Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0)[END_REF], represented in the form (46) along with the functions/matrices defined in (49) satisfy the constraints in [START_REF] Smyshlyaev | Adaptive boundary control for unstable parabolic PDEs -Part II: Estimation-based designs[END_REF]. The example also highlights that all terms (49) present in the condition are not zero. As a result of the satisfaction of the conditions of Theorem 2, we conclude that the PDE ( 43) is exponentially stable. We have thus shown the exponential stability for a backstepping controlled system in its original coordinates, highlighting the potential application of the proposed conditions to other boundary controlled PDEs for which a simple stable system in some target coordinates is not available.

V. CONVEX RELAXATION: PDES WITH POLYNOMIAL DATA

In the last section we formulated the stability conditions as the search for variables satisfying a non-linear PDME, subject to an ODMI and two LMIs. Even though we showed, via the backstepping boundary controlled PDE in [START_REF] Yamashita | Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0)[END_REF], that a choice for such variables which satisfy the conditions of Theorem 2 exists, finding such variables can be difficult in general. Therefore, it is of interest to formulate a computationally tractable numerical test to verify the conditions of Theorem 2. In this section we provide such a numerical test for PDEs of the form ( 8), but defined by polynomial data. We formulate a convex feasibility problem to verify the conditions of Theorem 2, and hence, establish the exponential stability of the zero solution of [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF].

Consider [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF] defined by the polynomial data

A1 ∈ R β×β(α+1) [x], A2, A3 ∈ R β×β(α+1) [(x, y)], (65a)
F1 ∈ R βα×β(α+1) [x], F2 ∈ R βα×2βα . (65b)
For such PDEs, we present the following corollary to Theorem 2.

Corollary 1. For any positive scalars ǫ, δ, polynomial degree d ∈ N, polynomial matrix Z q(α,β,d) ∈ R q(α,β,d) [(x, y)] defined in (1) and with p = q(α, β, d), r = q(0, β, d),

suppose there exist

T ∈ Π R 2 (r, 0, β), R ∈ Π R 2 (p, α, β), (66a)
{K1, K2, H1, H2, B1, B2} ∈ Π R 3 (α, β), (66b)
ST ∈ S β+2r [x], SR ∈ S β(3α+1)+2p [x], (66c)
S ∈ S β(3α+1) [x], (66d)
such that

T (x) - ǫI β 0 β,2r 0 2r,β 02r -ST (x)ω(x) ∈ Σ β+2r [x], (67a)
ST ∈ Σ β+2r [x], (67b)
R(x) -SR(x)ω(x) ∈ Σ β(3α+1)+2p [x], (67c)
SR(x) ∈ Σ β(3α+1)+2p [x], (67d) He (U b (x) + K b (x) + H b (x) + B b (x)) -R b (x) -S(x)ω(x) ∈ Σ β(3α+1) [x], (67e) S(x) ∈ Σ β(3α+1) [x], (67f) Ū (x, y) + H(x, y) + B(x, y) -R(x, y) = 0α+1, (67g)
where ω(x) = x(1 -x) and

{U b , Ū } =Ξ1(A1, A2, A3, T b , T , δ), {T b , T } =Ξ2(Z q(0,β,d) , T ), {R b , R} =Ξ2(Z q(α,β,d) , R), {K b , H b , H, B b , B} =Ξ3(K1, K2, H1, H2, B1, B2),
and where Πi, Ξi, i ∈ {1, 2, 3}, are defined in (2)- [START_REF] Bucci | Optimal boundary control with critical penalization for a PDE model of fluid-solid interactions[END_REF].

Then, the PDE (8) defined with polynomial data (65) is exponentially stable.

Proof: We begin by recalling that, as defined in Section I-C, the set of Sum-of-Squares (SOS) polynomial matrices Σ n [x] contains symmetric polynomial matrices S n [x] which are positive semidefinite for all x ∈ R. Furthermore, ω(x) = x(1 -x) satisfies ω(x) ≥ 0, for all x ∈ [0, 1]. Therefore, if (67a)-(67b) are satisfied, then

T (x) - ǫI β 0 β,2q(0,β,d) 0 2q(0,β,d),β 0 2q(0,β,d) 0, ∀x ∈ [0, 1].
Therefore, for the integral expression T (w) in ( 27) defined using (66a) and Vr(x, y) = Z q(0,β,d) (x, y), using Proposition 1 we conclude that there exist positive scalars θ1, θ2 such that

θ1 u 2 L 2 ≤ T (u) ≤ θ2 u 2 L 2 , ∀u ∈ L2 [0, 1]; R β . (68)
Similarly, if (67c)-(67d) hold, then the integral expression R (w) defined in ( 24) using (66a) and Yp(x, y) = Z q(α,β,d) (x, y), using Lemma 2 satisfies

R (u) ≥ 0, ∀u ∈ H α [0, 1]; R β . (69)
Additionally, (67e)-(67f) imply

He (U b (x) + K b (x) + H b (x) + B b (x)) -R b (x) 0, ∀x ∈ [0, 1]. (70)
Then, using (67g) and ( 68)-(70) we may follow the same line of reasoning as in the proof of Theorem 2 to conclude the assertion of this corollary. Unlike the conditions (36) in Theorem 2, the conditions in (67) are linear in the unknown variables (66) since we have fixed Yp(x, y) = Z q(α,β,d) (x, y) and Vr(x, y) = Z q(0,β,d) (x, y) to define R and T (see [START_REF] Kreyszig | Introductory functional analysis with applications[END_REF] and [START_REF] Krstic | Adaptive boundary control for unstable parabolic PDEs -Part I: Lyapunov design[END_REF], respectively). Moreover, the set of polynomials is closed under the operations of differentiation and integration. Therefore, since the variables in (66) are polynomial matrices, the expressions in (67) are polynomials. Thus, in order to analyze stability of (8) defined by polynomial data in (65), Corollary 1 requires the search of polynomial matrices which, under linear operations, need to belong to either the set of SOS polynomials as in (67a)-(67f), or satisfy affine constraints as in (67g). The search for SOS polynomials under affine constraints is a Semi-Definite Programming (SDP) feasibility problem [5, Chapter 3], [START_REF] Vandenberghe | Semidefinite programming[END_REF]. Since an SDP is a convex optimization problem, in order to establish the stability of (8) defined by (65), we can thus test the feasibility of the following convex optimization problem:

Find (66) subject to (67). (71)
To solve this problem, one may use freely available packages SOS-TOOLS [START_REF] Valmorbida | SOSTOOLS: Sum of squares optimization toolbox for MATLAB[END_REF] or YALMIP [START_REF] Lofberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF]. These packages simplify the extraction of the underlying SDP to (71) by providing an interface for the declaration of polynomial variables (66) and constraints (67). The SDP problem associated with (71) can then be solved by SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones[END_REF] or SDPA [START_REF] Yamashita | Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0)[END_REF].

A. Numerical Examples

We now test the presented numerical methodology on examples of PDEs of the form [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF] and defined by the polynomial data (65). In particular, we consider the examples provided in ( 9), ( 10), ( 12), [START_REF] Dullerud | A course in robust control theory: A convex approach[END_REF], and [START_REF] El-Farra | Analysis and control of parabolic PDE systems with input constraints[END_REF]. The studied examples depend on a scalar parameter λ and the stability of trajectories is guaranteed for λ ∈ [0, λ stable ). As stated earlier, such a parameterization allows us to verify the effectiveness of the proposed method by comparing the maximum value of λ for which (71) is feasible to the value of λ stable .

We perform the numerical experiments with ǫ = 10 -3 and δ = 10 -4 . Furthermore, polynomial degrees are kept to a maximum of 4 as this is a limitation imposed by the available memory of the computer on which these experiments are performed (a random access memory of 8 gigabytes). In order to search for the maximum λ for which (71) is feasible, we perform a bisection search with a resolution of 10 -3 .

Example 1: Partial (Integro) Differential Equation (P(I)DE): We begin by considering the P(I)DE provided in [START_REF] Coron | Dissipative boundary conditions for one-dimensional nonlinear hyperbolic systems[END_REF]. We use finitedifferences with spatial discretization of 1500 uniformly spaced points to approximate that ( 9) is stable for λ < 3.728. The P(I)DE ( 9) may be cast as [START_REF] Castillo | Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control[END_REF] with

α = 1, β = 1, w(x, t) = v(x, t), (72a) A1(x) = 0 1 , A2(x, y) = λ(x -y) 0 , (72b) A3(x, y) = λ(x + y) 0 , F1(x) = 0 βα,β(α+1) , (72c) F2 = 1 0 . ( 72d
)
Example 2: Wave Equation: Now consider the PDE in [START_REF] Coron | Stabilization of a rotating body beam without damping[END_REF]. From the solution [START_REF] Curtain | An introduction to infinite-dimensional linear systems theory[END_REF], it is evident that we must have λ ∈ (0, 1) for the exponential stability of [START_REF] Coron | Stabilization of a rotating body beam without damping[END_REF]. In fact, the system is finite-time stable for λ = 0.

To perform the stability analysis using the proposed method, let us cast [START_REF] Coron | Stabilization of a rotating body beam without damping[END_REF] as ( 8) by defining w(x, t) = w1(x, t) w2(x, t)

⊤ = ∂xv(x, t) ∂tv(x, t) ⊤ , yielding ∂tw(x, t) = ∂t ∂xv(x, t) ∂tv(x, t) = ∂t∂xv(x, t) ∂ 2 x v(x, t)
= ∂x w2(x, t) w1(x, t) , (73a) w1(0, t) -(1 -λ)w2(0, t) = 0, w2(1, t) = 0.

(73b)

The form (73) is also adopted in the analysis performed in [28, Section 7.1] and [START_REF] Fridman | An LMI approach to H∞ boundary control of semilinear parabolic and hyperbolic systems[END_REF]Section 3.1]. The main motivation for using (73) lies in the fact that the L2-norm of w is representative of the kinetic and potential energies of the trajectories of ( 10) by involving the L2-norm of ∂tv and ∂xv, respectively. In turn, this form can be cast as ( 8) by choosing

α = 1, β = 2, w(x, t) = w1(x, t) w2(x, t) ⊤ , (74a)
A1(x) = 0 0 0 1 0 0 1 0 , A2(x, y) = A3(x, y) = 0 β,β(α+1) , (74b)
F1(x) = 0 βα,β(α+1) , F2 = 0 0 1 -(1 -λ) 0 1 0 0 . ( 74c
)
Example 3: Parabolic PDE with Scalar Coefficients: Now consider the parabolic PDE in [START_REF] Datko | Extending a theorem of A. M. Liapunov to Hilbert[END_REF]. Using the solution of this PDE in [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled first-order hyperbolic linear PDEs with a single boundary input[END_REF], it is evident that this PDE is exponentially stable for λ < π 2 /4. The PDE (12) may be written in the form of (8) by setting where ξ and ν are real valued functions. With this representation the PDE ( 14) can be written as

α = 2, β = 1, w(x, t) = v(x, t), (75a) A1(x) = λ 0 1 , A2(x, y) = A3(x, y) = 0 β,β(α+1) , (75b) F1(x) = 0 βα,β(α+1) , F2 = 0 0 1 0 0 1 0 0 . ( 75c
∂tξ(x, t) = ∂ 2 x ν(x, t), ∂tν(x, t) = -∂ 2 x ξ(x, t), ∂xξ(0, t) + (1 -λ)ν(0, t) = 0, ∂xν(0, t) + (λ -1)ξ(0, t) = 0, ξ(1, t) = 0, ν(1, t) = 0,
which is a system of PDEs linearly coupled both in the domain and on the boundaries. With this representation, we obtain (8) by choosing

α = 2, β = 2, w(x, t) = ξ(x, t) ν(x, t) ⊤ , (76a)
A1(x) = 0 0 0 0 0 1 0 0 0 0 -1 0 , ( 76b
) A2(x, y) = A3(x, y) = 0 β,β(α+1) , F1(x) = 0 βα,β(α+1) , (76c)
F2 =     0 0 0 0 0 1 -λ 1 0 0 0 0 0 λ -1 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0     . (76d)
Example 5: Parabolic PDE with Distributed Coefficients and Boundary Feedback: Finally, let us consider the PDE in [START_REF] El-Farra | Analysis and control of parabolic PDE systems with input constraints[END_REF]. Using a finite-difference scheme with 1500 uniformly spaced spatial points we approximate that this PDE is exponentially stable for λ < 5.217.

We can cast [START_REF] El-Farra | Analysis and control of parabolic PDE systems with input constraints[END_REF] in the form of (8) by setting

α = 2, β = 1, w(x, t) = v(x, t), (77a) A1(x) = θ(x) ∂xφ(x) φ(x) , (77b) A2(x, y) = A3(x, y) = 0 β,β(α+1) , (77c)
F1(x) = ψ(x) 0 0 ξ(x) 0 0 , F2 = 0 0 0 1 1 1 0 0 . ( 77d
)
For each of the example PDEs in ( 9), ( 10), ( 12), [START_REF] Dullerud | A course in robust control theory: A convex approach[END_REF], and (15) represented in the form of ( 8) in ( 72), ( 74), ( 75), (76), and (77), respectively, we perform a bisection search on the parameter λ solving problem (71) for each fixed value of the parameter. The optimal obtained values λmax are presented in Table I. In each of the examples considered, the proposed methodology was able to prove the stability within 99.9% of calculated/approximated stability margin λ.

Other examples such as the coupled hyperbolic PDEs of the form considered in [START_REF] Di Meglio | Stabilization of a system of n + 1 coupled first-order hyperbolic linear PDEs with a single boundary input[END_REF] and the Euler-Bernoulli beam model may be found in [START_REF] Gahlawat | A semi-definite programming approach to stability analysis of linear partial differential equations[END_REF] wherein a preliminary result of the one in this work was presented.

VI. CONCLUSION AND FUTURE WORK

We presented a method to assess exponential stability of a large class of linear PDEs. The method is based on a Lyapunov Functional (LF) approach that reduces the stability analysis to the verification of integral inequalities on the subspaces defined by the boundary conditions of the PDEs. The verification of these inequalities is performed by solving a system of coupled differential equations and inequalities which contains a non-linear Partial Differential Matrix Equation (PDME), a linear Ordinary Differential Matrix Inequality (ODMI) and two Linear Matrix Inequalities (LMIs). The key element in relating the integral inequalities' verification to the solution of the system of coupled differential equations is the application of the Fundamental Theorem of Calculus (FTC) and Green's theorem.

We apply the proposed method to a backstepping controlled parabolic PDE from the literature. Using the solution to the system of coupled differential equations we obtain the LF certificates of stability.

Importantly, the proposed method is amenable to numerical formulations. Indeed, for the case of the class of PDEs defined by polynomial data, we provide a formulation in terms of a Sumof-Squares (SOS) program. The solution to the underlying convex optimization feasibility problem is obtained for examples of PDEs spanning the considered class of systems. We would like to highlight that the proposed methodology has been shown to work well on the considered numerical examples. Since the proposed conditions are only sufficient, a principled and exhaustive study should be performed to claim the effectiveness to all systems in the considered class.

The manuscript provides a proof of concept of our choice of LFs and the verification of the resulting integral inequalities. Therefore, we are working on the extension of the proposed framework to the problem of controller synthesis. In particular, for the large class of PDEs considered, we are interested in formulating methodologies for boundary and/or distributed controller synthesis. Further possible directions of this work include observer design using boundary/indomain state measurements, optimal control design, and stabilization of uncertain and non-linear systems.

+ ˆ1

x T (y, x) ⊤ w(y)dy.

Since the operator Ξ is self-adjoint, we may use (8) to obtain

∂tV (w) = 1 2 Ξ∂tw, w L 2 + 1 2 Ξw, ∂tw L 2 = 1 2 Ξw, ∂tw L 2 + 1 2 Ξw, ∂tw L 2 = Ξw, ∂tw L 2 = ˆ1 0 T b (x)w(x)dx + ˆx 0 T (x, y)w(y)dy + ˆ1 x T (y, x) ⊤ w(y)dy ⊤ × A1(x)wα(x)
+ ˆx 0 A2(x, y)wα(y)dy + ˆ1 x A3(x, y)wα(y)dy dx.

Therefore,

V d (w) = - 4 i=1 Φi, (85)
where, The term Φ1 may be written as

Φ1 = ˆ1 0 T b (x)w(x) + ˆx 0 T (x, y)w(y)dy + ˆ1 x T (y, x) ⊤ w(y)dy ⊤ A1(x)wα(x)dx, Φ2 = ˆ1 0 (T b (x)w(x)) ⊤ × ˆx 0 A2(x,
Φ1 = ˆ1 0 wα(x) ⊤ T b (x)A1(x) 0 β,2βα 0 3βα,β(α+1) 0 3βα,2βα wα(x)dx + ˆ1 0 ˆx 0 wα(x) ⊤ T (x, y) ⊤ A1(x) 0 βα,β (α+1)
⊤ wα(y)dydx

+ ˆ1 0 ˆ1 x wα(x) ⊤ T (y, x) ⊤ A1(x) 0 βα,β (α+1)
⊤ wα(y)dydx.

Then, applying Lemma A.1 to the double integrals and writing the single integral kernel in a symmetric form produces

Φ1 = ˆ1 0 wα(x) ⊤ He T b (x)A1(x) 0 β,2βα 0 3βα,β(α+1) 0 3βα,2βα wα(x)dx + ˆΩ wα(x) ⊤ Γ 1 2 T (x, y) ⊤ A1(x) 0 βα,β(α+1) ⊤ + 1 2 T (x, y)A1(y) 0 βα,β(α+1) wα(y)dydx. (86)
The term Φ2 may be written as

Φ2 = ˆ1 0 ˆx 0 wα(x) ⊤ T b (x)A2(x, y) 0 βα,β(α+1) wα(y)dydx + ˆ1 0 ˆ1 x wα(x) ⊤ T b (x)A3(x, y) 0 βα,β(α+1)
wα(y)dydx.

Applying Lemma A.1 produces

Φ2 = ˆΩ wα(x) ⊤ Γ 1 2 T b (x)A2(x, y) 0 βα,β(α+1) + 1 2 T b (y)A3(y, x) 0 βα,β(α+1) ⊤ wα(y)dydx. (87)
The term Φ3 may be written as

Φ3 = ˆ1 0 ˆx 0 T (x, y) ⊤ 0 βα,β ⊤ wα(y)dy+ ˆ1 x T (y, x) 0 βα,β ⊤ wα(y)dy × ˆx 0 A2(x, y)wα(y)dy + ˆ1 x A3(x, y)wα(y)dy dx. (88)
Then, applying Lemma A.2 to (88) with

F1(x, y) = T (x, y) ⊤ 0 βα,β ⊤ , F2(x, y) = T (y, x) 0 βα,β ⊤ , G1(x, y) =A2(x, y), G2(x, y) = A3(x, y), v(y) = wα(y), produces Φ3 = 1 2 ˆΩ wα(x) ⊤ Γ [U ] wα(y)dydx, (89)
where U (x, y) is defined in [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF].

Finally, the term Φ4 may be written as Substituting (86), (87), (89) and (90) into (85) produces [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF].

Φ4 = ˆ1 0 wα(x) ⊤ He δT b (x) 0 β,3βα 0
We now provide the proofs of the claims in the paper. We start by providing a proof of Lemma 1.

Proof of Lemma 1: Let us choose the LF candidate as V (w). Then, as shown in the beginning of this appendix, along the solutions of the system -∂tV (w(•, t)) -2δV (w(•, t)) = V d (w(•, t)).

Since (18b) holds, we have that for all solutions w of (8)

V d (w(•, t)) = -∂tV (w(•, t)) -2δV (w(•, t)) ≥ 0, ∀t ≥ 0.

Integrating this expression in time produces

V (w(•, t)) ≤ e -2δt V (w(•, 0)), and thus, using (18a) produces µ1 w(•, t) 2 L 2 ≤ e -2δt µ2 w(•, 0) 2 L 2 .

Then, we conclude that [START_REF] Gahlawat | A convex sum-of-squares approach to analysis, state feedback and output feedback control of parabolic PDEs[END_REF] Since R(x) 0, for all x ∈ [0, 1], we conclude that

R (u) = ˆ1 0 f (x) ⊤ R(x)f (x)dx ≥ 0, ∀u ∈ H α [0, 1], R β .
Proof of Proposition 1: It is easily established that the functions T b and T are continuous on their respective bounded domains of definitions. Therefore, there exists a scalar 0 < θ2 < ∞ such that ˆ1 0 u(x) ⊤ T b (x)u(x)dx + ˆΩ u(x) ⊤ Γ T (x, y) u(y)dydx

≤ θ2 u L 2 , (91)
for all u ∈ L2 [0, 1]; R β . Now, let us suppose that T (x) satisfies (28a), then using the proof of Lemma 2, it can be established that we get that ˆ1 0 u(x) ⊤ T b (x)u(x)dx + ˆΩ u(x) ⊤ Γ T (x, y) u(y)dydx

≥ ǫ u 2 L 2 , ∀u ∈ L2 [0, 1]; R β . (92)
Setting ǫ = θ1 and using (91) and (92) completes the proof for the case when the matrix T (x) satisfies (28a).

Let us now assume that T (x) satisfies (28b). Define the operators M and K as (Mu) (x) = Q1(x)u(x), (Ku) (x) = ˆx 0 Q2Vr(x, y)u(y)dy.

Since Q1(x) is continuous on [0, 1] and satisfies (28b), there exists a scalar µ > 0 yielding

ǫ 2 u 2 L 2 ≤ Mu 2 L 2 ≤ µ u 2 L 2 , ∀u ∈ L2 [0, 1]; R β . ( 93
)
The inverse of operator M is well defined since Q1(x) is invertible on [0, 1]. We thus have

M -1 Ku (x) = ˆ1 0 G(x, y)u(y)dy, (94)
where G(x, y) = Q1(x) -1 Q2Vr(x, y), x ≥ y 0, y > x .

From (28b) we have ˆ1 0 u(x) ⊤ T b (x)u(x)dx + ˆΩ u(x) ⊤ Γ T (x, y) u(y)dydx

= ˆ1 0 f (x) ⊤ T (x)f (x)dx = ˆ1 0 f (x) ⊤   Q1(x) ⊤ Q ⊤ 2 0 r,β     Q1(x) ⊤ Q ⊤ 2 0 r,β   ⊤ f (x)dx, = (M + K) u, (M + K) u L 2 = M I + M -1 K u, M I + M -1 K u L 2 , ( 95
)
where I is the identity operator. Owing to the continuity of Q1(x) -1 and Vr(x, y), it is straightforward to verify that

ˆ1 0 ˆ1 0 |G(x, y)| dydx < ∞.
Therefore, the operator M -1 K in (94) is compact on L2 [0, 1]; R β [START_REF] Renardy | An introduction to partial differential equations[END_REF]Theorem 7.83]. Since M -1 K is compact, from [16, Theorem 5, Appendix D] we have that the range of the operator I + M -1 K is closed. Moreover, it can be established using [START_REF] Kreyszig | Introductory functional analysis with applications[END_REF] that the null space of the operator I + M -1 K is the set {0} ∈ L2 [0, 1]; R β . Since the operator I + M -1 K is linear and its nullspace is {0} ∈ L2 [0, 1]; R β , we have that it is injective. Since I + M -1 K has a closed range and is injective, using [2, Theorem 2.5] we conclude that there exists a scalar γ > 0 such that

(I + M -1 K)u 2 L 2 ≥ γ u 2 L 2 , ∀u ∈ L2 [0, 1]; R β . ( 96
)
Therefore, using (96) and ( 93) we get

M I + M -1 K u, M I + M -1 K u L = M I + M -1 K u 2 L 2 ≥ ǫ 2 I + M -1 K u 2 L 2 ≥ ǫ 2 γ u 2 L 2 , ∀u ∈ L2 [0, 1]; R β .
Substituting this expression into (95) and setting θ1 = ǫ 2 γ completes the proof. The proof for the case when T (x) satisfies (28c) follows similarly.

Proof of Lemma 4: Consider the vector field φ1(x, y) φ2(x, y) = uα-1(x) ⊤ H1(x, y)uα-1(y) uα-1(x) ⊤ H2(x, y)uα-1(y) .

Then, by Green's theorem ˛∂Ω (φ1(x, y)dx + φ2(x, y)dy)

+ ˆΩ (∂yφ1(x, y) -∂xφ2(x, y)) dydx = 0, where ∂Ω denotes the boundary of the domain Ω. Then, the proof is completed by using the definition of the vector field, the definitions of the projection matrices in Section I-C, and by applying Lemmas A.1 and A.2.

Proof of

Example 4 :

 4 Complex Valued In-Domain and Boundary Coupled Parabolic PDE: Consider the PDE ∂tv

  and j is the imaginary unit. The well-posedness of this PDE is established using [25, Theorem 3.1]. Example 5: Parabolic PDE with Distributed Coefficients and Boundary Feedback: Finally, let us consider the following PDE with boundary feedback and polynomial coefficients ∂tv

  y)wα(y)dy + ˆ1 x A3(x, y)wα(y)dy dx, Φ3 = ˆ1 0 ˆx 0 T (x, y)w(y)dy + ˆ1 x T (y, x) ⊤ w(y)dy ⊤ × ˆx 0 A2(x, y)wα(y)dy + ˆ1 x A3(x, y)wα(y)dy dx, Φ4 =δ ˆ1 0 w(x) ⊤ T b (x)w(x)dx + δ ˆΩ w(x) ⊤ Γ T (x, y) w(y)dydx.

  ˆΩ wα(x) ⊤ Γ 2δ T (x, y) 0 β,βα 0 βα,β 0 βα wα(y)dydx. (90)

  holds with κ = µ2/µ1. Proof of Lemma 2: If we define f (x) =   ūα(x) ´x 0 Yp(x, y)uα(y)dy ´1 x Yp(x, y)uα(y)dy   , then, using Lemma A.1 and Lemma A.2, it can be shown that R (u) = ˆ1 0 f (x) ⊤ R(x)f (x)dx.

-ǫ u 2 L 2 =

 22 ) ⊤ T b (x)u(x)dx + ˆΩ u(x) ⊤ Γ T (x, y) u(y)dydx ˆ1 0 f (x) ⊤ T (x) -´x 0 Vr(x, y)u(y)dy ´1 x Vr(x, y)u(y)dy   . Now, since T (x) -ǫI β 0 β,2r 0 2r,β 02r 0, ∀x ∈ [0, 1],

Lemma 5 :

 5 Since for all u ∈ B,ˆ1 0 F1(x) F2 ūα(x)dx = 0 βα,1 ,we get that for all u ∈ B ūα(x)dx= ˆ1 0 ˆ1 0 ūα(y) ⊤ B1(y) B2 F1(x) F2 ūα(x)dydx = ˆ1 0 ˆ1 0 ūα(y) ⊤ B1(y)F1(x) B1(y)F2 B2F1(x) B2F2 ūα(x)dydx.Then, using the fact that ūα(x) = uα(x) u b α ⊤ and applying Lemma A.1 completes the proof.

)

  Example 4: Complex Valued In-Domain and Boundary Coupled Parabolic PDE: Now consider the PDE in (14). This PDE is exponentially stable for λ < 1 [28, Exercise 6.2, Eqns. 6.74-6.76].

		λ stable	λmax	% Accuracy
	Example 1 (9)	3.728	3.725	99.91
	Example 2 (10)	1	0.999	99.90
	Example 3 (12)	π 2 /4	2.466	99.94
	Example 4 (14)	1	0.999	99.90
	Example 5 (15)	5.217	5.215	99.96
	We may write

v(x, t) = ξ(x, t) + jν(x, t),

TABLE I :

 I Maximum λ, λmax, for which the problem (71) is feasible for the example PDEs. Here, λ stable is the analytically/numerically determined stability margin and percentage accuracy is calculated as (λmax/λ stable ) × 100.

For brevity we have dropped the temporal dependency of w.

VII. ACKNOWLEDGMENTS

The authors would like to express their gratitude to Dr. Federico B. Argomedo for providing assistance with the numerical simulation of Partial (Integro)-Differential Equations P(I)DEs.

APPENDIX A

In this Appendix we provide a few results which are used in the manuscript.

Lemma A.1. For any α, β ∈ N and Lebesgue integrable functions K1 : Ω → R β(α+1)×β(α+1) and K2 : Ω → R β(α+1)×β(α+1) , the following identity holds ˆ1 0 ˆx 0 uα(x) ⊤ K1(x, y)uα(y)dydx

The proof is established in a straightforward manner by applying Fubini's theorem (change of order of integration) followed by a switch between variables x and y.

Lemma A.2. For any u ∈ L2 [0, 1]; R β , α, β ∈ N, and Lebesguemeasurable and L ∞ functions F1, G1 : Ω → R β×β , F2, G2 : Ω → R β×β , the following identity holds

where

Proof: We begin by observing that the

where the linear bounded operators on L2 [0, 1]; R β are defined as

and where the Hilbert adjoint of the operator F is given by

Therefore, we get

For each of the double integral terms defining (F ⋆ Gv) (x), we change the order of integration and switch between the variables y and z and substitute the result into (79). The proof is then completed by applying Lemma A.1.

Remark. Lemma A.2 also holds for any

and with u(x) replaced by uα(x) in (78).

Lemma A.3. For the function J in (44), h1 in (50), h2 in (51), b in (59), and g in (62), the following equation holds for any scalar δ

where we have used the fact that

Now, we have that J satisfies (45a), i.e.,

Therefore, we also have the following

We have from (45b) that

and we also have that

Combining the above expressions gives

Substituting ( 82) into (81) produces

Applying the FTC (integration by parts) to the last term gives Substituting into (83) produces

APPENDIX B

In this Appendix we provide the proofs of previously stated results. We begin by showing that for the solution w of the PDE (8)

V d (w(•, t)) = -∂tV (w(•, t)) -2δV (w(•, t)), where V (u) and V d (u) are defined in ( 16) and [START_REF] Fantuzzi | Semidefinite relaxation of a class of quadratic integral inequalities[END_REF], respectively. Let us write the integral expression in (16) as 1

where the self-adjoint operator Ξ on L2 [0, 1]; R β is defined as (Ξw) (x) =T b (x)w(x) + ˆx 0 T (x, y)w(y)dy