OBLIQUE PROJECTION BASED STABILIZING FEEDBACK FOR NONAUTONOMOUS COUPLED PARABOLIC-ODE SYSTEMS

. Global feedback stabilizability results are derived for nonautonomous coupled systems arising from the linearization around a given time- dependent trajectory of FitzHugh–Nagumo type systems. The feedback is explicit and is based on suitable oblique (nonorthogonal) projections in Hilbert spaces. The actuators are, typically, a ﬁnite number of indicator functions and act only in the parabolic equation. Subsequently, local feedback stabiliz- ability to time-dependent trajectories results are derived for nonlinear coupled parabolic- ode systems of the FitzHugh–Nagumo type. Simulations are presented showing the stabilizing performance of the feedback control.


1.
Introduction. In this work we focus on coupled systems consisting of a parabolic equation with controls, and an ordinary differential equation. This class of systems include well-known models describing electric excitations and the propagation of electric waves in nerve fibers and in heart tissue. In this context the system of equations is referred to as the monodomain equations. Given a bounded domain Ω ⊂ R d , where d is a positive integer, with smooth boundary Γ = ∂Ω, the controlled monodomain equations read with Neumann boundary conditions Finally u = u(t), taking values in R M , is a control function at our disposal.
In electrophysiology, see [11,Section 12.3.3], the variable v models the transmembrane electric potential of the human heart and w represents a gating variable. Typical models include the FitzHugh-Nagumo, the Rogers-McCulloch, and the Aliev-Panfilov model. See [7,18,21,2,8]. The coupling M takes the form where the constants d, e, γ, ρ are also given, and are all strictly positive with the following exceptions defining the model: (2b) Assume that a desired heart rhythm is given as the solution of the (uncontrolled) system (1), corresponding to a suitable initial condition (v 0 , w 0 ) ∈ L 2 (Ω) × L 2 (Ω), If (v 0 , w 0 ) = (v 0 , w 0 ), then the behavior of (1) without control can be considerably different from the desired behavior of (v, w), even if (v 0 , w 0 ) − (v 0 , w 0 ) is small. Here we look for a feedback control u, so that the solution of (1) goes to the desired heart rhythm (v, w), solving (3), provided the difference (v 0 , w 0 ) − (v 0 , w 0 ) is small.
We will follow a standard idea: first we find a feedback operator which stabilizes globally the linearization of (1) around (v, w), then we use a suitable fixed point argument to conclude that the same feedback operator also stabilizes the nonlinear system locally.
More precisely, by direct computations, we can see that the (controlled) linearization around (v, w), satisfies ∂ ∂t y + Ay + A r y + Sz − Without loss of generality we may suppose that the set of actuators {1 ωi | i ∈ {1, 2, . . . , M }} ⊂ L 2 (Ω) is linearly independent. We also set U M := span{1 ωi | i ∈ {1, 2, . . . , M }}. Notice the subscript r in the reaction operator A r , in order to distinguish it from the diffusion operator A.
For simplicity let us denote L 2 := L 2 (Ω) and H 1 := H 1 (Ω). We will show that a sufficient condition for global stabilizability of system (4), and for local stabilizability of system (5), is given by α M +1 > 6 + 4 P Contents. The rest of the paper is organized as follows. In Sect. 2 we derive results on the stability of linear coupled systems in an abstract setting. These results are applied in Sect. 3 to obtain stabilizability conditions for general linear coupled parabolic-ode systems. In Sect. 4 we derive the stabilizability result for the concrete example of the linearized monodomain equations, in particular, it contains the proof of Main Theorem 1.1 in Sect. 4.1 and that of Main Theorem 1.2 in Sect. 4.2. Numerical simulations are presented in Section 5, for both the linearized system (9) and the nonlinear system (10), confirming the theoretical results and showing the performance of the explicit oblique projection stabilizing feedback control. Finally the Appendix gathers comments on the existence and uniqueness of weak solutions for the systems involved in the main text.
Notation. We write R and N for the sets of real numbers and nonnegative integers, respectively, and we define R r := (r, +∞), for r ∈ R, and N 0 := N \ {0}.
We denote by Ω ⊂ R d a bounded open connected subset, with d ∈ N 0 . For a normed space X, we denote by | · | X the corresponding norm, by X its dual, and by ·, · X ,X the duality between X and X. The dual space is endowed with the usual dual norm: |f | X := sup{ f, x X ,X | x ∈ X and |x| X = 1}. In case X is a Hilbert space we denote the inner product by (·, ·) X .
For an open interval I ⊆ R and two Banach spaces X, Y , we write W (I, X, Y ) := {f ∈ L 2 (I, X) | ∂ ∂t f ∈ L 2 (I, Y )}, where the derivative ∂ ∂t f is taken in the sense of distributions. This space is endowed with the natural norm |f | W (I, X, Y ) := 1/2 . In case X = Y , we write H 1 (I, X) := W (I, X, X).
The time derivative (in the distribution sense) of a vector function v taking values in a Banach space X will be denoted byv := d dt v. If the inclusions X ⊆ Z and Y ⊆ Z are continuous, where Z is a Hausdorff topological space, then we can define the Banach spaces X × Y , X ∩ Y , and X + Y , endowed with the norms |(a, b)| X×Y := |a| 2 X + |b| 2 Y 1 2 ; |a| X∩Y := |(a, a)| X×Y ; and |a| X+Y := inf (a X , a Y )∈X×Y |(a X , a Y )| X×Y | a = a X + a Y , respectively. We can show that, if X and Y are endowed with a scalar product, then also X × Y , X ∩ Y , and X + Y are. In case we know that X ∩ Y = {0}, we say that X + Y is a direct sum and we write X ⊕ Y instead.
Again, if X and Y are endowed with a scalar product, then also W (I, X, Y ) is. The space of continuous linear mappings from X into Y will be denoted by L(X, Y ). When X = Y we simply write L(X) := L(X, X).
if the inclusion is also dense, respectively compact. C [a1,...,a k ] denotes a nonnegative function of nonnegative variables a j that increases in each of its arguments.
2. Stability of coupled evolutionary systems. We shall present a stability result for a coupled abstract system in a form which is convenient for our purposes.
We shall utilize a separable Hilbert space H, which we consider as a pivot space, H = H . We introduce a closed subspace G ⊆ H, and two evolutionary systemṡ in the Hilbert spaces G and H, where for suitable Hilbert spaces V ⊂ G and W ⊂ H, we suppose that A(t) ∈ L(V, V ) and D(t) ∈ L(W, W ) are operators in G and H, respectively, with domains D(A) = D(A(t)) ⊆ G and D(D) = D(D(t)) ⊆ H, independent of t. Furthermore, we assume that Since the existence and uniqueness of solutions is not the focus of this paper, it is simply assumed by us. In the Appendix, we briefly recall the standard procedure on the construction of weak solutions as a weak limit of suitable Galerkin approximations. In particular we will see that the above assumption is satisfied for a general class of operators A, D, S, and R. Applications to the linearized monodomain equations will be given in Sect. 4.
holds true, then for all ε sufficiently small there exists a constant D c ≥ 1 such that the weak solution of system (15) satisfies independently of (v 0 , w 0 ).
The proof is given below. Before we derive some auxiliary results, where Assumptions 2.2 and 2.3 are assumed to hold true.
Proof. Using Duhamel formula, we integrate the equation in (15b) and obtain, for t ≥ s 0 , Therefore, for t ≥ s 0 , we obtain Let us consider first the case µ A = µ D . In this case we obtain Now, from (23), (24), and (25), it follows It remains to prove that the last inequality also holds in the case µ A = µ D . Thus, let µ A = µ D and notice that we have U −D (t,s0) w 0 any ρ ∈ (0, µ D ), t ≥ s 0 ≥ 0. We know, from (26), that from which we conclude that (26) also holds in the case µ A = µ D .
Corollary 2.6. If inequality (18) holds true, then the weak solution of system (15) satisfies for a suitable constant D c ≥ 1, independent of (v 0 , w 0 ).
We also have that z solvesż = −Az, z(s 0 ) = z 0 , if, and only if,

Now if (18) holds true, then it also holds
for a small enough ε ∈ (0, min{µ A , µ D }). From Corollary 2.6 it follows that for a suitable constant D c ≥ 1 independent of (v 0 , w 0 ), from which we can conclude (19).
3. Stability of the coupled parabolic-ode closed-loop system. Here we show that the explicit oblique projections based feedback control proposed in [13] for stabilization of nonautonomous linear parabolic equations is also able to stabilize a general class of nonautonomous linear coupled parabolic-ode systems, where the control acts (only) in the parabolic component. We will prove that the systeṁ with the explicit feedback is stable, under suitable assumptions on the linear span of the actuators U M = span{1 ω1 , 1 ω2 , . . . , 1 ω M } and on the operators in (33).
3.1. Assumptions. We look at system (33) as an evolutionary system in H × H, where H is our pivot separable Hilbert space, H = H. First we present our assumptions on the operators in (33). They will guarantee that the assumptions in the abstract setting of Section 2 are fulfilled. Most of the assumptions are standard, and concern the existence and uniqueness of weak solutions for auxiliary systems. The only nonstandard assumption is the main stabilizability condition given below in Assumption 3.7. It can be seen as a generalization of the condition presented in [13] for the case of parabolic equations.
Let V be another Hilbert space with V ⊂ H.
From now we will suppose that V is endowed with the scalar product (y, z) V := Ay, z V , V , which still makes V a Hilbert space. Necessarily, A : V → V is an isometry.
Assumption 3.2. The inclusion V ⊆ H is dense, continuous, and compact.
Necessarily, we have that y, z V , V = (y, z) H , for all (y, z) ∈ H × V, and also that the operator A is densely defined in H, with domain D(A) satisfying Further, A has a compact inverse A −1 : H → D(A), and we can find a nondecreasing system of (repeated) eigenvalues (α n ) n∈N0 and a corresponding complete basis of eigenfunctions (e n ) n∈N0 : We can define, for every β ∈ R, the fractional powers A β , of A, by For the time-dependent operators we assume the following: (14). For any s 0 ≥ 0 and w 0 ∈ H there is one, and only one, (H, W)-weak solution foṙ and there are constants C D ≥ 1 and µ D > 0, independent of (s 0 , w 0 ), such that the solution w satisfies |w(t) Note that in case K U M (t) = 0, we have that (33) is in the form of system (15), with G = H and V = V . Note also that (34) is in the form of system (15) We suppose that E ⊥ M ⊂ H and E M ⊂ H are endowed with the norm inherited from H, and we denote (cf. Lemma 2.5) and R := |R| L ∞ (R0,L(H)) .
Assumption 3.7. The linear span U M of the actuators satisfies , (36a) .

3.2.
The stability result. The main result of this section, which will lead to main Theorems 1.1 and 1.2, is the following.
Theorem 3.10. Under Assumptions 3.1-3.7, the coupled system (33) is stable: there are constants C c ≥ 1 and µ c > 0 such that The proof is presented below, in section 3.2.2.

Auxiliary results.
Here we derive some auxiliary results we will use in the proof of Theorem 3.10. We suppose that the Assumptions 3.1-3.7 do hold true. Notice that G = E ⊥ M is a closed subspace of H. Let us set V = G ∩ V , which we endow with the norm inherited from V . Observe that V is a closed subspace of V , and A maps V onto V .
Moreover, we have Proof. The existence and uniqueness is proven in [13, Section 3.1]. Note that we and for any given positive constants γ 1 and γ 2 , and from (35b), it follows that we can choose γ 1 and γ 2 such that which implies (38).
If there is a pair of constants C q ≥ 1 and λ q > 0 such that Furthermore, the constants D q and ε q are independent of (s 0 , v 0 , w 0 ). The constants µ and µ D are as in (38) and Assumption 3.4.
Proof. When q = 0, the existence and uniqueness of the solution, in the space where µ is as in (38). Then, we observe that for small enough ε ∈ (0, min{µ A , µ D }) we have Observe also that system (40) is in the form of system (31), by setting η = Rq Then from Corollary 2.7, it follows (41), for all ε q ∈ (0, min{ε, λ q }) and a suitable constant D q ≥ 1.

4.
Applications. Stabilization of the monodomain equations. Assume that the desired trajectory (v, w) is given as the solution of the uncontrolled system (1). As mentioned in the Introduction, the difference (y, z) := (v, w)−(v, w) satisfies system (10), with the operators defined as in (6). In this section we prove Theorems 1.1 and 1.2 which concern the stabilization of systems (9) and (10)   The proof follows from standard arguments. A few details are given in the Appendix. Thus, under condition (8) and Assumption 3.7, all the assumptions in Section 3.1 are satisfied. By Theorem 3.10 this implies the following Corollary, which in turn implies Theorem 1.1. From Remark 3.8 we recall that condition (7), which is assumed in Theorem 1.1, implies that Assumption 3.7 is satisfied.
Corollary 4.2. Let λ > 0 and s 0 ≥ 0. If (8) and Assumption 3.7 hold true, then there are constants C ≥ 1 and µ > 0, such that the solution of the linear system (9) satisfies Here the constants C and µ are independent of the triple (s 0 , y 0 , z 0 ).

4.2.
Proof of main Theorem 1.2. To deal with the nonlinear systems, we will need strong solutions. We can prove that such solutions exist due to the fact that the reaction operator A r defined as in (6) satisfies if (8) The proof, which relies in part on known arguments, is sketched in the Appendix. The next lemma gathers estimates on our nonlinearity N = N (y) = N (t, y). Let us write, for simplicity, Further for any given ς > 0 there exists C 2 ≥ 0 such that with {ε 1 , ε 2 } ∈ [0, +∞) and {ε 3 , ε 4 , ε 5 , ε 6 } ∈ [2, +∞).
Proof. Recall that from (6) Concerning G 4 we proceed as follows: estimate (43) is trivially satisfied for e = 0, thus we consider only the case e = 0, then with (ζ y , ζ z ) := p −p = (y −ỹ, z −z) ∈ D(A) × H, we find Next we consider the case d = 3, that is, Ω ⊂ R 3 . An analogous argument can be followed for d ∈ {1, 2}. To further estimate these two inequalities, we use the for any given ς > 0, and for suitable constants C 1 > 0, C 2 > 0 and C 3 = C 3 (ς) > 0. Therefore we can conclude that the inequalities in (43) also hold true for G 4 . Finally, concerning G 3 we proceed as follows: First H . Thus (43) holds true also for N replaced by G 3 . and from (8) and (6) we have A r ∈ L ∞ (R 0 , L(H, V )) and S ∈ L ∞ (R 0 , L(H)).
Remark 4.6. We are now prepared to provide the proof for Theorem 1.2, which relies on fixed point arguments. With two modifications, we can rely on the arguments utilized in [5, Section 3.2] (see also [19,Section 3]). First, the feedback operator in these references is obtained as the solution of a suitable differential Riccati equation, in place of K U M as in (33c). However, thanks to the estimates in Section 3] are slightly stronger than those available here.
In fact, in [19,Section 3], instead of property (43b), the following stronger assumption is assumed: Since Assumption 3.7 is implied by (7), Theorem 1.2 follows from the following.
Proof. We follow the main argument sketched in [19,Section 3] and [3, Section 4]. Let us consider the systeṁ  L(H,H)) . Here F (1,1) y := K U M (y, 0) and F (1,2) z := K U M (0, z). We will look for a solution of the nonlinear system in a subset Z µ ⊂ Z µ of the Banach space

KARL KUNISCH AND SÉRGIO S. RODRIGUES
For a given constant > 0 we define the subset Z µ as follows.
Further we define the mapping Ψ : Z µ → Z µ loc ,v → v, taking a given vectorv to the solution v ofv s Step [ii]: Existence. The fixed point argument. Proceed as in [19,Section 3], by using (50), with N (v) in the place of f , and using (43a), we can conclude that Therefore, we can conclude that if z 0 ∈ V is sufficiently small, |z 0 | 2 V < , then there exists a unique fixed point z = Ψ(z) =z ∈ Z for Ψ. That is, Necessarily y z := v solves (10). Further, notice that (45) follows from z ∈ Z µ ρ . s Step [iii]: Uniqueness. We show now that a solution for (10) in the space Z := L 2 loc (R 0 , D(∆)) C([0, +∞), V) ⊃ Z µ is unique. Indeed, let z 1 and z 2 be two solutions in Z, for (51). Then e := z 1 − z 2 solves (49) with g = N (z 1 ) − N (z 2 ) in the place of N (z). Using (43b), and following standard arguments, we can obtain d dt |e| 2 Using e(s 0 ) = 0 we find z 2 (t) − z 1 (t) = e(t) = 0, for t ≥ s 0 . This ends the proof.

Numerical simulations.
Here we present simulation results for (9) and (10). As actuators we will take suitable (piecewise constant) indicator functions {1 ωi | i ∈ {1, 2, . . . , M }} ⊂ H. In this section H := L 2 (Ω). In Figure 1 we can see the placement of (rectangular) actuators' regions. Recall also that the possibly repeated eigenvalues α i of the Laplacian, under Neumann boundary conditions increase likely linearly, as we would expect from Weyl asymptotic formula [24]. See also [12,10].
The following simulations were carried out using a finite elements based spatial approximation of our equations, with a Crank-Nicolson time discretization.
We will set the parameters in the FitzHugh-Nagumo, Rogers-McCulloch, and Aliev-Panfilov models so that we will have 3 constant steady states ( v i , w i ), i ∈ {1, 2, 3} and so that the first (pde) components are v 1 = 0, v 2 = 1, and v 3 = 2. As we will see:   ( We will argue that for all 3 models: • the linearizations around the steady state ( v 1 , w 1 ) and around the steady state ( v 3 , w 3 ) are stable, • the linearization around the steady state ( v 2 , w 2 ) is unstable. We will consider the systems (9) and (10) with two choices as reference trajectories. The first one is the unstable steady state the second one is spatially and temporally dependent and is defined as follows: we consider the time-dependent vector function T = (T 1 , T 2 ) : R 0 → R 2 , where s ∈ N stands for the smallest integer which is smaller than s ≤ s < s + 1. Then we define our reference trajectory (v(t, x), w(t, x)) as with T 1 (t, x 1 , x 2 ) := 0.001 sin(πt) cos(πt + 0.1( Notice that (v, w) is periodic in time with period 2, (v(t + 2, x), w(t + 2, x)) = (v(t, x), w(t, x)) for all t ≥ 0. Furthermore, for any given nonnegative integer m ∈ N, (v(2m, x), w(2m, x)) = ( v 2 , w 2 ) and (v(2m + 1, x), w(2m + 1, x)) = ( v 3 , w 3 ). That is, the trajectory (v, w) "oscillates" between the unstable steady state ( v 2 , w 2 ) and the stable steady state ( v 3 , w 3 ).
The initial condition in all the simulations, unless otherwise explicitly stated, is the constant In the feedback (33c), we have set the parameter 5.1. The linearized FitzHugh-Nagumo model. Here we consider the linear system (9), with the parameters in (6) set as We can see that the vectors as in (52a) are constant steady states of the uncontrolled and unforced nonlinear coupled system (1) (i.e., with u = 0 and f = 0). Indeed we can see that for a constant vector ( v, w) we have ∆ v = 0 and that Therefore with the parameters as in (56) we find that ξ(0, 0) = ξ(1, 10) = ξ(2, 20) = 0 0 , which shows that the vectors as in (52a) are constant steady states for system (1). Observe also that the linearization (4) is stable around (0, 0) and (2,20), and is unstable around (1,10). Indeed, the linearization around a constant vector ( v, w) reads and we find that with the parameters as in (56): . It is clear that both eigenvalues have strictly negative real part, because (1 + δ) 2 − 4(δ + dγ) < (1 + δ) 2 .
In the figures below the annotation "dyn=feed" refers to the dynamics under the feedback control and "dyn=free" to the free dynamics (without control). The annotation "model=..." specifies the type of the monodomain model defined by the selection of parameters (cf. (2)).
In Figure 2 we see the performance of our feedback control, which is able to stabilize system (9). The slope of log(|(y, z)| Recall that the explicit feedback has been constructed to stabilize the pde component, and as a corollary the stability of the ode component follows as well. Consequently we cannot expect to obtain a rate of decay for the ode component which is better (smaller) than −δ.
In Figure 2 we can see also the curve illustrating the squared norm Further notice that the free dynamics is exponentially unstable, the norm is increasing exponentially with rate approximately 1 2 120 30 = 2.

5.1.2.
The case of a time-dependent trajectory. We consider the time-dependent trajectory (v, w) as in (53b). In Figure 3 we observe that after time t ≥ 5 the norm of the solution of the linear system (9) decreases exponentially with rate approximately −δ = −0.01. The free dynamics is exponentially unstable with norm increasing exponentially with rate approximately 1.

5.2.
The nonlinear FitzHugh-Nagumo model. We consider again the timedependent reference trajectory as in (53b), but now for the nonlinear system (10).

5.3.1.
The case of the constant unstable steady state. Here the reference trajectory (v, w) = (1, 10) is as in (53a). In Figure 5 we can see that the explicit feedback is able to stabilize the system (9) exponentially with rate approximately −δ = −0.01, while the free dynamics is exponentially unstable with rate approximately 2.  Figure 5. Norm of solution. Linearization around the steady state.

5.3.2.
The case of a time-dependent trajectory. We consider the time-dependent trajectory (v, w) as in (53b). In Figure 6 we can see that the explicit feedback is able to stabilize the system (9) exponentially with rate approximately −δ = −0.01, while the free dynamics is exponentially unstable with rate approximately 1.5.

The nonlinear Rogers-McCulloch model.
Here, we consider the timedependent trajectory (v, w) as in (53b), and consider the corresponding nonlinear system (10) with the parameters as in (59). In Figure 7 we see that our explicit feedback is able to exponentially stabilize system (10). The norm of the controlled system decreases exponentially with rate approximately −δ = −0.01. Notice also that the free dynamics is unstable.
Note that from (57), with the parameters as in (60), we find which shows that the vectors as in (52b) are constant steady states for system (1). For the linearization around a constant steady state (58), we find that with the parameters as in (59): • For ( v, w) = (0, 0), the eigenvalues of the matrix −c −d γ −δ are {−c, −δ}.
They are both strictly negative.
5.5.1. The case of the constant unstable steady state. Here we consider the timeindependent trajectory (v(t), w(t)) = (1, 5), is an in (53a). In Figure 8 we see that our explicit feedback is able to stabilize system (9). The norm of the solution decreases exponentially with rate approximately −δ = −0.01. The norm of the free dynamics is exponentially increasing with rate approximately 0.5.

5.5.2.
The case of a time-dependent trajectory. We consider again the linear system (9), with the time-dependent trajectory (v, w) as in (53b). In Figure 9 we confirm that also in this case the feedback is able to stabilize system (9). The norm decreases exponentially with rate approximately −δ = −0.01. The free dynamics is not stable. The norm of the solution is not converging to zero, and is "close" to a periodic behavior.  5.6. The nonlinear Aliev-Panfilov model. We take the time-dependent trajectory (v, w) as in (53b), and the corresponding nonlinear system (10), again with the parameters as in (60). For Figure 10 we ran the simulation for time t ∈ [0, 30] as in previous examples. In this situation the results do not allow us to conclude that the feedback is stabilizing. The norm seems to decrease but, we cannot yet clearly see a rate of stability. Therefore we ran the simulation for a larger time interval and the results are shown in Figure 11. Again we observe that an exponentially decreasing rate is achieved. In fact, it is approximately −δ = −0.01 (the slope for time t ≥ 100 is approximately −0.02).
We also observe that the norm corresponding to the free dynamics is not converging to zero, that is the free dynamics is not stable.

5.7.
On the sign of the leading coefficient a. In the above examples the linearization based feedback control is able to stabilize the nonlinear system for the initial condition (−1, −1), as in (54). Notice that in the models above, the sign of the leading coefficient a = 1 of av 3 − b 2 v 2 + cv is positive, and in this case we may say that is has the good sign, in the sense that such sign somehow favors stability of the nonlinear system, for example it may guarantee the existence of global solutions (cf. [9, Chapter 5, Remark 5.1], for the uncoupled parabolic equation). The positive sign of a may explain why the feedback still works for the nonlinear system for such a "large" initial condition.
Recall that simulations in [19], for a single uncoupled parabolic equation, show that the feedback is able to stabilize the nonlinear system only if the initial condition is small enough. This is also the statement of Theorem 4.7, concerning the coupled system.
Below we will consider the "a−" model where we take a negative leading coefficient a. Notice that our results are independent of the sign of a. We will consider the following parameters which roughly can be seen a variant of the Aliev-Panfilov model with the negative sign for a and c and with b = 0. In such situation we may expect the uncontrolled solution of the uncoupled corresponding parabolic system to blow up in finite time. See [9,Chapter 5], [14], and references therein. Therefore we may expect the solution of the corresponding coupled system to explode as well. 5.7.1. The linearized "a−" model. We consider the time-dependent trajectory (v, w) as in (53b), with (1,1) in the place of the steady states: that is with ( v 2 , w 2 ) := (1, 1) =: ( v 3 , w 3 ).
In Figure 12 we see that the norm of the controlled solution of system (9) decreases exponentially with rate approximately −δ = −0.01. The norm of the free dynamics increases exponentially with rate approximately 4. 5.7.2. The nonlinear "a−" model. Here we consider the nonlinear system (10). We will test with the resized initial conditions First in Figure 13, for the free dynamics for the nonlinear system (10), We observe that the norm of the solution blows up in finite time, for all initial conditions as in (62). Next we consider the corresponding systems under the action of the feedback control. In this case, in Figure 14, we can see that our feedback is not able to stabilize the system for the initial conditions (−1, −1) and (−0.1, −0.1). But, it is able to stabilize the system if we take the smaller initial condition (−0.05, −0.05) which is consistent with our local stabilization result for the nonlinear system, as stated in Theorem 1.2. 6. Final remarks.
6.1. On the stability conditions (7). We have shown that stability of the explicit closed-loop system holds provided the conditions in (7)  Condition (7b), on the other hand is not so transparent. But for rectangular domains we have investigated whether, for large enough M , the condition can be satisfied for M actuators 1 ω M i , i ∈ {1, 2, . . . , M }, whose supports are constrained to cover altogether an a priori fixed arbitrarily small volume, vol with r ∈ (0, 1) independent of M . This question was answered affirmatively for 1D, in [20,Theorems 4.4,4.5,and 5.2], for both Dirichlet and Neumann boundary conditions. Therefore, from [13, Section 4.8.1], we also know that we can place suitable M d actuators in higher-dimensional rectangles in R d , d ≥ 2, so that an analogous condition to (7b) is satisfied, which again implies the stability of the closed-loop system.
Finally for general (nonrectangular) domains Ω ⊂ R d , d ≥ 2, we do not know whether we can construct M actuators 1 ω M i (covering an a priori fixed total volume) so that (7b) is satisfied. For this, it would be enough to prove that it is possible to construct/place such actuators in such a manner that the operator norm of the oblique projection P E ⊥ M U M will remain bounded as M increases, because α M +1 → +∞. This possibility was conjectured in [13,Section 4.8.2]. This is an interesting open problem, because simulations in [13] and in this manuscript show the stabilizing performance of the feedback in nonrectangular 2D domains. 6.2. On robustness. One of the reasons for the importance of feedback control in applications is its gain in robustness against noise and disturbances. Simulations in [13], for a single parabolic equation, show that our feedback control can attenuate state estimation errors of the form η(t, x) = η 1 (t, x)y(t, x) + η 2 (t, x), where η 1 (t, x) and η 2 (t, x) reach small magnitudes. In particular, the component η 1 (t, x)y(t, x) can be seen as a disturbance of the reaction operator A r . In addition to the numerical evidence, it would be interesting to investigate the robustness properties of the explicit proposed feedback, due to the importance of such properties for automatic control applications. Also other types of noise and plant (dynamics) disturbances could be of interest, because the type of such errors may depend on the problem/applications. 6.3. On the number of actuators. Here we consider the same setting as in Section 5.7.1. Recall that in Figure 12, we have the results corresponding to the case of 28 actuators as in Figure 1. Now in Figure 15, we present the corresponding results in the case we take only 4 actuators as in Figure 15. We constructed the actuators so that the volume covered by the 4 actuators in Figure 15 is equal to the total volume covered by the 28 actuators in Figure 1. In either case the total volume covered by the actuators is approximately 5.4% of the volume of Ω. Though the volume covered by the actuators is the same, we conclude that the explicit feedback corresponding to the 4 actuators as in Figure 15 is not able to stabilize system (9), while the explicit feedback corresponding to the 28 actuators as in Figure 1 is, as shown by Figure 12. From the results in [13] we can conclude that once we know that P remains bounded with respect to M , then increasing both λ and the number of actuators allows to achieve any prescribed rate of exponential stability in the case of a single uncoupled parabolic equation. In Figure 16 we observe that 4 actuators chosen as in Figure 15 are not able to stabilize the uncoupled parabolic equation, while 28 actuators as in Figure 1 are able to stabilize the uncoupled parabolic equation with rate approximately 1, which is the best rate we can get, independently of M , because we have chosen λ = 1 in our feedback law, see (55).   As we said before, besides the number of actuators, another important question concerns the (optimal) placement of the actuators, say so that P is as small as possible. In Figure 1 and Figure 15 we simply tried to "spread" the actuators over Ω as "uniformly" as possible. However, the results in [20], show that the best way of "spreading" the actuators is a nontrivial question already in the 1D case.
By curiosity, the norm (computed numerically, see [13]) of the projection with the 4 actuators as in Figure 15 is given by P Even if two different placements of the actuators lead to the stability of the closed-loop linear system, a better placement of the actuators may lead to a better performance of the feedback control, see [13, Section 5.1].
6.4. On the stabilization to periodic solutions. While our results apply for general bounded trajectories, in our numerical examples we have taken periodic trajectories as targeted solutions. We will show here that for this case our feedback law is also periodic. Indeed, suppose that our targeted solution (v, w) for system (3) is time-periodic: for a given period T > 0 (v(s + T ), w(s + T )) = (v(s), w(s)), for all s ≥ 0. (63) In that case we observe that the feedback law is also periodic Indeed, for all y, z ∈ H × H, Ay + A r (s)y − λy + S(s)z = K U M (s)(y, z), because, from (6), we see that A r (s+T )y = A r (s)y and S(s+T ) = S(s), when (v, w) is T -periodic, as in (63).
Recalling that the monodomain equations provide a model for the heart rhythm we expect the healthy solution to be periodic, that is, it makes sense to consider at the very beginning that (v, w) satisfies (63).
Then it is natural to ask whether we can take advantage of the time periodicity of A + A r + S.
In [4,17] the stabilization of time-periodic parabolic systems is proven under suitable conditions on the set of actuators. Such conditions take advantage of the fact that the operators in the parabolic equation are time-periodic, which allows to derive suitable properties of the asymptotic behavior of the solutions from the spectral properties of the periodic map or Floquet map.
The feedback law proposed in [4,17] is, however, not explicit. It is given implicitly and involves the solution of a suitable periodic Riccati equation. The computation of such a solution can be an expensive numerical task.
It would be interesting to know whether we can construct an ad-hoc, still explicit, feedback operator in this situation. Recall that in [13, Section 6.5], an ad-hoc explicit feedback is proposed for the case of time-independent targeted trajectories (i.e., steady states).
which we can obtain after multiplying (68) by z N , imply that with D 2 independent of N and z 0 . Hence, by (68), we obtain ż N 2 L 2 ((s0,s0+t),V ) = Q N Lz N 2 ,s0+t),V) , with D 3 independent of (N, s 0 ). Here we use that the orthogonal projection Q N ∈ L(H) is also defined and orthogonal in V (cf. [13,Lemma 3.3]), which implies that Q N 2 L(V ) = 1. Note that for two eigenpairs (α i , e i ), (α j , e j ) of A = −ν∆ + Id we have (e i , e j ) V := (A − 1 2 e i , A − 1 2 e j ) L 2 = (A −1 e i , e j ) L 2 = α −1 i (e i , e j ) L 2 . Now, it is straightforward to check that points (2) and (3), in Definition A.2, hold true.
Since the assumption in Lemma A.3 is satisfied, it follows the existence and uniqueness of a weak solution for system (65), with L as in (67). This finishes the proof of Proposition 4.1.