Kernel Well-Posedness and Computation by Power Series in Backstepping Output Feedback for Radially-Dependent Reaction-Diffusion PDEs on Multidimensional Balls

Recently, the problem of boundary stabilization and estimation for unstable linear constant-coefﬁcient reaction-diffusion equation on n -balls (in particular, disks and spheres) has been solved by means of the backstepping method. However, the extension of this result to spatially-varying coefﬁcients is far from trivial. Some early success has been achieved under simplifying conditions, such as radially-varying reaction coefﬁcients under revolution symmetry, on a disk or a sphere. These particular cases notwithstanding, the problem remains open. The main issue is that the equations become singular in the radius; when applying the backstepping method, the same type of singularity appears in the kernel equations. Traditionally, well-posedness of these equations has been proved by transforming them into integral equations and then applying the method of successive approximations. In this case, with the resulting integral equation becoming singular, successive approximations do not easily apply. This paper takes a different route and directly addresses the kernel equations via a power series approach (in the spirit of the method of Frobenius for ordinary differential equations), ﬁnding in the process the required conditions for the radially-varying reaction (namely, analyticity and evenness) and showing the existence and convergence of the series solution. This approach provides a direct numerical method that can be readily applied, despite singularities, to both control and observer boundary design problems.

In the early and mid-1990s, inspired by "bifurcation control," wildly popular in the physics community at that time, and linked to interesting applications, Krener and Kang developed control designs for normal forms that include locally quadratic dependencies, which remain the definitive results of the art on this subject.
Art Krener cast an eye on PDE control systems at least as early as the mid-1990s, in the framework of a US Air Force funded project on nonlinear control of jet engine instabilities. While other researchers focused on the first-order Galerkin approximations of the models of rotating stall instabilities in axial flow compressors in jet engines, Art formulated and studied generalized higherorder Moore-Greitzer nonlinear models [15].
It is a delight to see Art Krener take on PDE control as the preoccupation for the present stage of his career [20,21]. Going beyond the conventional development of operator Riccati equations, and the limitation to the study of their well posedness, in his trademark fashion, Art is producing computable approximate solutions to Riccati PDEs using Al'brekht's approach, which he has already brought to its state-of-the art form for nonlinear and stochastic ODEs.

Introduction
In this paper we introduce an explicit boundary output-feedback control law to stabilize an unstable linear radially-dependent reaction-diffusion equation on an n-ball (which in 2-D is a disk and in 3-D a sphere).
This paper extends the spherical harmonics [7] approach of [36], which assumed constant coefficients, using some of the ideas of [40]. For a finite number of harmonics, we design boundary feedback laws and output injection gains using the backstepping method [23] (with kernels computed using a power series approach) which allows us to obtain exponential stability of the origin in the L 2 norm. Higher harmonics will be naturally open-loop stable. The required conditions for the radially-varying coefficients are found in the analysis of the numerical method and are non-obvious (evenness of the reaction coefficient). The idea of using a power series to compute backstepping kernels was first seen in [4] (without much analysis of the method itself, but rather numerically optimizing the approximation) and later in [10], where piecewise-smooth kernels require the use of several series. Here, we prove that the method provides a unique converging solution, in the spirit of the method of Frobenius for ordinary differential equations.
Some partial results towards the solution of this problem were obtained in [38] and [37] for the disk and sphere, respectively; however they required simmetry conditions. Older results in this spirit were obtained in [34] and [28]. This paper extends and completes our conference contribution [41] where the ideas where initially presented (without proof).
To the best of our knowledge, this paper presents the first rigorous proof of convergence of a power series solution for the backstepping kernel equations. Thus, this work consolidates the method as a valid alternative to more traditional numerical approaches, which include finite difference approximations of the kernel equations [23,17,11,3], the use of symbolical successive approximation series [35], or the numerical solution of the integral version of the kernel equations [16,6]. The main advantages of the method are its simplicity (it does not require the sometimes cumbersome conversion to integral equations thus preventing mistakes or any consideration about discrete meshes), speed (modern computing systems can reach high orders of the series in seconds), precision (one reaches a simple polynomial in one variable for the gain at the boundary that does not require interpolation), adaptability (it can be adapted to settings with discontinuous kernels by breaking the domain in pieces, see [10]), and capacity to produce kernels depending on parameters (by symbolically solving the kernel equations). The main drawback is the analyticity requirement of the system coefficients, even though most physical systems and examples seen in backstepping papers indeed posses analytic coefficients, and possibly a slow convergence rate of the series in some cases.
The backstepping method has proved itself to be an ubiquitous method for PDE control, with many other applications including, among others, flow control [32,39], nonlinear PDEs [33], hyperbolic 1-D systems [14,5,2], or delays [24]. Nevertheless, other design methods are also applicable to the geometry considered in this paper (see for instance [31] or [8]).
The structure of the paper is as follows. In Section 3 we introduce the problem. In Section 4 we state our stability result. We study the well-posedness of the kernels in Section 5, which is the main result of the paper, proving existence of the kernels and providing means for their computation; interestingly, odd and even dimensions require a slightly different approach. We briefly talk next about the observer in Section 6, but skip most details based on its duality with respect to the controller. Then, we give some simulation results in Section 7. We finally conclude the paper with some remarks in Section 8.

n-D Reaction-Diffusion System on an n-ball
Consider the following constant-coefficient reaction-diffusion system in an n-dimensional ball of radius R: where u = u(t, x), with x = [x 1 , x 2 , . . . , x n ] T , is the state variable, evolving for t > 0 in the n-ball B n (R) defined as with boundary conditions on the boundary of B n (R), which is the (n − 1)-sphere S n−1 (R) defined as The boundary condition is assumed to be of Dirichlet type, where U(t, x) is the actuation. On the other hand the measurement y(t, x) is defined as where ∂ r denotes the derivative in the radial direction (normal to the (n − 1)-sphere), which would be defined as ∂ r u(t, x) = ∇u · x x . Differently from [36], we consider non-constant λ( x) verifying the following assumption.
is an analytic function of x and depends exclusively on the radius r = x .
Following [36], both the state and the actuation variable can be written in n-dimensional spherical coordinates, also known as ultraspherical coordinates (see [7], p. 93), which consist of one radial coordinate r and n − 1 angular coordinates θ. Then, using a (complex-valued) Fourier-Laplace series of Spherical Harmonics 1 to handle the angular dependencies, defined as where N(l, n) is the number of (linearly independent) n-dimensional spherical harmonics of degree l, given by N(0, n) = 1 (representing the mean value over the n-ball) and, for l > 0, with Y n lm being the m-th n-dimensional spherical harmonic of degree l. The coefficients in (6)-(7) are possibly complex-valued.
Following [36] and using (6)-(7) one reaches the following independent complex-valued 1-D reaction-diffusion equation for each spherical harmonic coefficient: evolving in r ∈ (0, R], t > 0, with boundary conditions In these equations, we have considered Dirichlet boundary conditions. The measurement would be the flux at the boundary, namely ∂ r u m l (t, R). Note that following [12, p. 640], a second boundary condition, reflecting the second-order character of (9) and the need to avoid singular behaviors, can be expressed as: In the above equations, the integers m and l stand for the order and degree of the harmonic, respectively. Note that the higher the degree (corresponding to high frequencies), the more "naturally" stable (9)-(10) is, as seen next. Define the L 2 norm and the associated L 2 space as usual, where | f | 2 = f f * , being f * the complex conjugate of f . Lemma 3.1. Given λ(r) and R, there exists L ∈ N such that, for all l > L, the equilibrium u m l ≡ 0 of system (9)- (10) is open loop exponentially stable, namely, for U m l = 0 in (10) there exists a positive constant D 1 , such that for all t u m l (t, ·) L 2 ≤ e −D 1 t u m l (0, ·) L 2 .
D 1 is independent of l, and only depends on n, λ(r), ε, and R, and can be chosen as large as desired just by increasing the values of L.
The proof is skipped as it mimics [40] just by using the L 2 norm as a Lyapunov function and Poincare's inequality.
Thus one only needs to stabilize the unstable mode with l < L. Since the different modes are not coupled, it allows us to stabilize them separately and re-assembling them. Moreover since only a finite number of harmonics is stabilized, there is no need to worry about the convergence of the control law as in [36], with its Spherical Harmonics series being just a finite sum.
Our objective can now be stated as follows. Considering only the unstable modes, design an output-feedback control law for U m l using, for each mode, only the measurement of ∂ r u m l (t, R). Our design procedure is established in the next section along with our main stability result.

Stability of controlled harmonics
Next, for the unstable modes we design the output-feedback law. The observer and controller are designed separately using the backstepping method, by following [36]; in this reference it is shown that both the feedback and the output injection gains can be found by solving a certain kernel PDE equation, which is essentially the same for both the controller and the observer. Thus, for the sake of brevity and to avoid repetitive material, we only show how to obtain the (full-state) control law, giving the basic observer design and some additional remarks later in Section 6.

Design of a full-state feedback control law for unstable modes
Based on the backstepping method [23], our idea is utilizing an invertible Volterra integral transformation where the kernel K n lm (r, ρ) is to be determined, which defined on the domain T k = {(r, ρ) ∈ R 2 ; 0 ≤ ρ ≤ r ≤ R} to convert the unstable system (9)-(10) into an exponentially target system: where the constant c > 0 is an adjustable convergence rate. From (14) and (16), let r = R, we obtain the boundary control as the following full-state law Following closely the steps of [36] to find conditions for the kernels, and defining K n lm (r, ρ) = G n lm (r, ρ)ρ ρ r l+n−2 , we finally reach a PDE that the G-kernels need to verify: with only one boundary condition: We assume as usual that these kernel equations are well-posed and the resulting kernel is bounded in T ; this will be analyzed later in Section 5, providing also a numerical method for its computation.

Closed-loop stability analysis of unstable modes
To obtain the stability result of closed-loop system, we need three elements. We begin by stating the stability result for the target system. We follow by obtaining the existence of an inverse transformation that allows us to recover our original variable from the transformed variable. Then we relate the L 2 norm with spherical harmonics. With these elements, we construct the proof of stability mapping the result for the target system to the original system. This is done by showing that the transformation is an invertible map from L 2 into L 2 .
We first discuss the stability of the target system, having the following lemma: where the constant D 2 is independent of n, l or m, and only depends on c, ε, and R; it can be chosen as large as desired just by increasing the value of c.
Proof. Consider the Lyapunov function: then, taking its time derivative, we obtaiṅ we then obtain, independent of the value of n, thus proving the result. For |l|≤ L, let c be chosen as in Lemma 4.1, and assume that the kernel K n lm (r, ρ) is bounded and integrable. The system (9) with boundary control (17) is closed-loop exponentially stable, namely there exists positive constants C and D 2 such that C and D 2 are independent of m or l, and only depend on n, L, λ(r), ε and R.
Proof. The proof consists of two parts, one is existence of an inverse transformation, and then showing the equivalence of norms of the variables u n lm and w n lm ; the result then follows from the stability of the target system.
As shown in [36], when K n (r, ρ) is bounded and integrable, the map (14) is invertible and its inverse transformation is which is also bounded and integrable. Call nowK andL the maximum of the bounds of the functioň K n lm andĽ n lm for a given n and all l ≤ L in their respective domains. It is easy to get where Combining then Lemma 4.1 with the norm equivalence between u m l and w m l system stated as in (27) and (28), it is easy to obtain Let C = √ M 1 M 2 , the result then follows.
Note that combining Lemmas 3.1 and 4.2 and taking D = min{D 1 , D 2 }, we get the following stability result for all spherical harmonics and thus the full physical system. Theorem 1. Under the assumption that the kernel K n lm (r, ρ) is bounded and integrable, the equilibrium u m l ≡ 0 of system (9)-(10) under control law (17) is closed-loop exponentially stable, namely, there exists a positive constant D, such that for all t where D can be chosen as large as desired just by increasing the value of L and c in the control design process.

Well-posedness of the kernel equations
Next, we state the main result of the paper, which was in part assumed in Theorem 1, also giving the requirements for λ(r). In addition the proof of the result also provides a numerical method to compute the kernels, which is an alternative to successive approximations which do not work in this case (due to the singularities at the origin; see for instance [38] to see the resulting singular integral equation that needs to be solved).

Theorem 2.
Under the assumption that λ(r) is an even real analytic function in [0, R], then for a given n > 1 and all values of l ∈ N, there is a unique power series solution G n lm (r, ρ) for (18)- (19), even in its two variables in the domain T , which is real analytic in the domain. In addition, if λ(r) is analytic, but not even, then there is no power series solution to (18)- (19) for most values l ∈ N.
The requirement of evenness for λ(r) might seem unusual. However, if we carefully consider Assumption 3.1, and since r = x = x 2 1 + x 2 2 + . . . + x 2 n , in physical space λ( x) will be nonanalytic, unless it is even. Thus, while solutions to the kernel equations might exist for non-even λ(r), we cannot expect them to be analytic. This result notwithstanding, if one is interested in controlling only very low-order harmonics, kernels do exist without this requirement, as shown in [38,37], which only consider the 0-th order harmonic (the mean) respectively for a disk and a sphere, and only require boundedness of λ(r).

Proof of Theorem 2
We start by giving out an algorithmic method to compute the power series for G n lm (r, ρ), which will allow us to prove Theorem 2 as well as numerically approximating the kernels.
First of all, we show that the evenness of λ(r) is a necessary condition to find an analytic solution. Next, it is possible to establish that the series for G n lm (r, ρ) only has even powers. Exploiting this property to suitably express (18)- (19), we finally show the existence of the power series and thus Theorem 2 follows. Convergence and related issues (radius of convergence) is studied towards the end, finishing the proof.

Computing a power series solution for the kernels
Starting from the most basic assumption of Theorem 2, we consider that λ(r) is analytic in [0, R], therefore it can be written as a convergent series (encompassing c and ε for notational convenience): which, by the evenness of λ, may only contain even powers 2 , this is, λ i = 0 if i is odd. We then, in the spirit of the method of Frobenius for ordinary differential equations, seek for a solution of (18)- (19) of the form: where the dependence on n, l and m has been omitted for simplicity (the solution will depend on these values). The series in (32) collects together (in the parenthesis) all the polynomial terms with the same degree. It is easy to see that the boundary condition (19) implies: which in particular implies C 00 = − λ 0 2ε . On the other hand, the left-hand side of (18) becomes where we have defined Finally, to express the right-hand side of (18), denote γ = n + 2l − 2 ≥ 0 and define the operators and thus, rewriting the sum to be homogeneous with (34),we find where, (assuming Equating (38) and (34), we obtain a system of equations: With λ(r) and n are fixed, we want to show that the kernel equations are solvable for all values of l ∈ N. Thus, γ takes increasing values. In addition we can assume γ = 1, since the case n = 3, l = 0 was already addressed in [37] showing that it reduces to the usual 1-D kernel equations for parabolic systems [23], which admits a power series solution according to [4].
The first two equalities, if γ = 1, imply and, in particular, C 10 = C 01 = 0, whereas the second equality results in a system of equations that needs to be solved recursively, starting at i = 0. It can be rewritten as follows to start at i = 2 (since C 00 , C 10 and C 01 are already determined).
Note that for each i ≥ 2, there are i + 1 coefficients in (32) but i + 2 relations: one from (33), two from (42) and i−1 from (43). Thus, it would seem that (33)-(42)-(43) is in general an incompatible system. This is indeed the case if λ(r) is not even, i.e., if the series (31) contains odd powers, as shown in the next section.

Evenness requirement of λ(r)
We start with the following result.
Lemma 5.1. If λ(r) is not even, then there are values of l ∈ N for which there is no solution to (18)- (19) in the form of (32).
Proof. We show that, if there exists i odd such that λ i = 0, then there is no solution in the form of a power series. First, if λ 1 = 0, then from (33) we know that C 01 +C 10 = − λ 1 4ε , however since form (42) one has C 01 = C 10 = 0, this cannot hold. Consider now there is indeed a value i > 1 for which a coefficient λ i is distinct from zero and let us show the result by contradiction. Consider the first such i. Now, since the right-hand side of (43) depends on C (i−2) j , one gets that for all odd i < i C i j must zero from (33)-(42)-(43) all having a zero right-hand side (this can be formalized with an induction argument; we skip the details). Thus, at i, the following system of equations has to be verified: and for 0 ≤ j ≤ i − 2, Let us consider l sufficiently large such that γ > i, so that the coefficient (46) is distinct from zero in the full range of j, namely 0 ≤ j ≤ i − 2. Then none of the coefficients in (46) is zero. Therefore, combining (44) with (46), from C i1 we can find C i3 , then C i5 , and so on. Similarly, from C i(i−1) we can find C i(i−3) , C i(i−5) and so on. These two sequences don't overlap because i is odd and therefore, one finds C i j = 0 for all 0 ≤ j ≤ i which is not compatible with (45) unless λ i = 0, which contradicts our initial assumption.
Next we show that evenness of λ implies evenness of the kernels.
Proof. We need to prove that C i j = 0 if either i or j is odd. From the proof of Lemma 5.1, we directly know that for odd i one has C i j = 0. Fix, then, i even and consider j odd; for i = 2, the result is obvious. Assuming C i j = 0 for all even numbers i < i and j odd, let us prove the result by induction on the first coefficient. As before, we would need to solve (45)-(43). The right-hand side (43) is zero as in (46) by the induction hypothesis (if k even) or directly zero if k odd. Then, following again the proof of Lemma 5.1, we have the same system of equations (45)-(46) for our even i and odd j's. Now: so starting from C i(i−1) = 0 we find C i(i−3) = 0, then C i(i−5) , and so on; however, with i being even, this sequence ends now in C i1 (thus, the proof of Lemma 5.1 does not apply because the sequences starting at C i1 and C i(i−1) overlap). Thus, one finds C i j = 0 for all odd values of j between 1 and i − 1.

Well-posedness of the coefficient system
Next, we show that the coefficients of the power series can always be found, which by the previous lemmas only requires studying the even coefficients. For simplification, we redefine (31) and (32) as: without bothering to redefine the coefficients (note that (35) does not require any change). Defining as well γ = γ 2 = n 2 + l − 1 ≥ 0, the new system of equations to be solved is and (49) Let us outline the solution procedure, and later derive some conclusions. Solving in (49) every C i j as a function of C i( j+1) we get: which can be written more briefly if we define, for i > 0 and 0 ≤ j < i, as To be able to simplify a bit the equation, redefinê then, and iterating this equality until reaching C ii , we get and inserting this into (48), we reach an equation for C ii , namely where and It is quite clear that these κ(i, γ ) will play an important role; in particular, if they are non-zero, one can always find a unique solution for the coefficients C i j . Thus one needs to show that κ(i, γ ) = 0 for any possible i or γ . The following lemma shows this is indeed the case, by exploiting a connection of the a i j coefficients with Gauss' hypergeometric functions.
Proof. Recalling from (56) and (51) the definitions of κ(i, γ ) and a i j , respectively, one has which can be rewritten in terms of binomial numbers and rising/falling factorials 3 [18] as and reordering the sum and using i j where the obvious fact (x) j = (−1) j (−x) j has been used. Consider now the finite polynomial p i (x, γ ) defined as From the definition of Gauss' hypergeometric function [1, p.561], denoted as 2 F 1 (a, b; c; x), in the polynomial case (a or b non-positive integer) and noting (−1) j i j = (−i) j j! , it is immediate that and therefore, from Gauss' summation theorem [1, p.556], which is applicable in this case since 1 + 2i > 0, finishing the proof.
The next result is an immediate conclusion of the positivity of κ(i, γ ): Lemma 5.4. If λ(r) is even, then, for all values of l ∈ R, the coefficientes in (47) that solve (18)- (19) can be uniquely found up to any order i.
To conclude the proof of Theorem 2, we need to prove analyticity of the series (47). This step, however, requires splitting the problem in two possible cases: odd dimension (thus, γ = n/2 + 2l − 1 is not an integer) and even dimension (γ integer).

Proof of analyticity for odd dimension
In the odd-dimension case, define the following coefficients: with L i0 defined as 1; these are well-defined given that γ is non-integer. They can also be expressed as Now, in (48)-(49), denote C i j = L i jČi j . Replacing in the recurrence we get Define nowB and the new set of recurrence equations forČ i j becomes rather simple: and the recurrence is easily solvable in terms of one element; for instance,Č ii : for j = 1, . . . , i. Replacing in (68) we reach Call Now, solving for the remaining coefficients from (70): Finally, recovering the coefficients C i j from C i j = L i jČi j and using (67): which is quite explicit. Notice that since ρ ≤ r, thus, defining α i = ∑ i j=0 |C i j |, if we can prove that ∑ ∞ i=0 α i r 2i converges for a certain radius of convergence R, so does G n lm (r, ρ) for ρ ≤ r ≤ R, and thus we obtain the required analyticity. Now: To prove the convergence of the power series ∑ ∞ i=0 α i r 2i consider the following lemma, inspired by [25]: Lemma 5.5. Consider g(x) = ∑ ∞ i=0 g i x 2i and h(x) = ∑ ∞ i=0 h i x 2i analytic functions, both with radius of convergence R. Let i 0 be a nonnegative integer, let (a i ) ∞ i=0 be a sequence of real numbers, and define f (x) = ∑ ∞ i=0 a i x 2i , where a i verify, for i > i 0 where the sequences b i , c i ≥ 0 are decreasing for i > i 0 , with c i also verifying lim i→∞ c i = 0. Then, f (x) is analytic with radius of convergence at least R.
Proof. Since g and h analytic with radius of convergence R we can write |g i |, |h i | ≤ MR −2i , where the definition as power series of squares has been taken into account. Thus: Defineǎ i = a i for i ≤ i 0 and, for i > i 0 ,ǎ i = b i MR −2i + c i ∑ i−1 j=0ǎ j MR −2i+1+2 j . Obvioulsy a i ≤ǎ i and therefore the radius of convergence of f (x) would be at least the radius of convergence of f (x) = ∑ ∞ i=0ǎ i x 2i . Now: where the inequality holds for sufficiently large i > i 0 and thus in the limit, therefore proving the lemma.
To apply Lemma 5.5 to (76) we need to bound some of the terms. In particular, if we are able to find b i and c i such that we get Thus, assuming that b i and c i verify the conditions given in Lemma 5.5, and given that λ(x) has a radius of convergence of at least R, we see that G n lm (r, ρ) converges and defines an analytic function for ρ ≤ r ≤ R, thus proving Theorem 2 for the odd-dimension case.
It remains to find such b i and c i . Proceeding exactly as in Lemma 5.3 with a slight modification, we directly find Now, let N = γ − 1/2 and i > 2N. One can see that for 1 ≤ j ≤ N, and Thus, for i > 2N, calling d i the following sequence it is clear that d i is a decreasing sequence, since from the ratio test r = lim i→∞ d i+1 It is obvious that b i is decreasing, since d i is decreasing. Now we need to find a sequence c i for the second term in (81). First of all, The following lemmas are needed to find a bound to (88).
Proof. Consider the ratio It is easy to see that Thus, the sequence always increases as long as j ≤ N, and we can look for a maximum j * > N. Then, for j > N, denote the ratio of (89) by f : Manipulating the expression, we find (i 2 − 1) − 2 j(i + 1) + γ (i + 1) ≤ 0 and canceling the term (i + 1) the following inequality for j is reached: Therefore, if (and only if) the bound given by (92) is verified, Therefore, we conclude that the maximum of the sequence |L i j | is reached at thus finishing the proof.
Thus, we are left with showing that |L i j * | |R i | is decreasing, which is expected, since L i j * is one of the elements that appear in the sum R i . Using the expression (83) and the formula for L i j that involves the Gamma function, we obtain the following: Now, the decreasing character of the sequence is established as follows (consider the case where i + N is odd, so that j * = i−1+N 2 ; the even case is analogous). Consider Stirling's approximation to the factorial, namely n! ≈ √ 2πn n e n . Then: Putting together (95)-(96), we obtain the following.
where we have broken the approximation into three functions: Notice that clearly lim i→∞ f 1 (i) = 0 (since f 1 (i) behaves like O(1/i) for large i), lim i→∞ f 2 (i) = 1, and lim i→∞ f 3 (i) = 16, thus it only remains to compute lim i→∞ f 2 (i) i , which is an indeterminate of the kind 1 ∞ . Resolving it (the details are omitted for brevity) one obtains that the limit is indeed 1. Thus, it is possible to find the decreasing sequence c i in (81), concluding the proof of convergence and analyticity in odd dimension.

Proof of analyticity for even dimension
The fact that γ is an integer makes the odd approach a priori impossible, since (65) would not be well defined (it would contain divisions by zero). However, to overcome that difficulty, we employ a partial solution for the kernel equations, to the order γ − 1, which helps to regularize the problem.
For this, consider F(r, ρ) = ∑ γ −1 i=0 r 2i φ i (ρ 2 ). Replacing this function in (18) and one gets the following recursive set of ODEs: which is solved starting at i = γ − 1: This can be written as which is an ODE with a regular singular point at x = 0. By applying the Frobenius method [13,Chapter 36] one can rewrite this equation as and its indicial equation is r(r − 1) + (1 + γ /2)r = 0, thus r 1 = 0 and r 2 = γ /2 (non-integer). We are interested in the solution of the form φ γ /2 = ∑ ∞ i=0 a i ρ 2i and discard the other solution. By Fuchs' theorem [9, p.146] this solution is analytic where λ(x) is analytic, thus the radius of convergence of the resulting φ γ −1 (ρ 2 ) is greater than one. Next, for i = γ − 2 up to i = 0: which, has the same indicial equation and again, also admits a solution in the required form. Applying once more Fuchs' theorem, this solution is analytic in intervals where both λ(x) and φ i+1 are analytic. Thus, by induction, we find a family of solutions such that the radius of convergence of all φ i is greater than R. The solutions just found have a degree of freedom (the first coefficient a i of their power series, which is φ i (0)). The idea is to construct the solution such that the boundary condition G n lm (r, r) = H(r) is satisfied up to order 2γ − 2. Thus: F(r, r) = ∑ γ −1 i=0 r 2i φ i (r 2 ) and expanding in power series φ i (r 2 ):
It can be shown that this scheme produces valid initial values for the φ i 's. However, an easier approach is to follow the general series approach of Section 5.1.1 up to order i = γ −1. By the uniqueness of the series development and identifying coefficients, it can easily be shown that φ i (0) = C ii . Next, calling G n lm (r, ρ) =Ǧ n lm (r, ρ) + F(r, ρ) the new boundary condition for the PDE becomes: G n lm (r, r) = H(r) − F(r, r) which starts at order 2γ . Thus, we can proposeǦ n lm (r, ρ) = r γ F 2 (r, ρ). One can see that the PDE for F 2 is and following previous sections, calling ψ(r 2 ) = H(r)−F(r,r) r 2γ and abusing the notation by keeping the same name for the coefficients C i j , one can find a power series development F 2 (r, ρ) = Now the approach of Section 5.1.4 becomes applicable and even easier, since all coefficients are positive. Indeed, define Mimicking Section 5.1.4 we reach where the last step can be performed due to the positivity of the redefined coefficients L i j compared to Section 5.1.4. Again, we apply Lemma 5.5 to (110). In this case, we define which is already a decreasing sequence. Then we need to find c i such that max r∈{0,...,i−1} and c i should be shown to be decreasing (for sufficiently large i) and convergent to zero. Consider the following lemma.
Lemma 5.8. Consider L i j as defined in (109) and R i , R i j as defined in (109). Then: the first property becomes evident, whereas the second property is immediate from the first since For the third property, note that The fourth property is obvious, noting that L i j ≤ L i( j+1) for j < i/2. Finally, for the last property, first note that it is only required to study 0 ≤ r ≤ i/2 given the third property. Now: and since this is an increasing function of r for 0 ≤ r < i, we can bound it by its value at r = i/2, thus proving the final property.
Therefore, setting c i = 4 i+2 i(i+2γ ) , a sequence that decreases to zero, we can apply Lemma 5.5 to (110) and follow the same steps as in Section 5.1.4 to obtain the result of Theorem 2 for the case of even dimensions. which can be written as 0 = λ(r) + ε∂ r P n lm (r, r) + ε d dr (P n lm (r, r)) + (n − 1) εP n lm (r, r) r + ε∂ ρ P n lm (r, r) − (n − 1) εP n lm (r, r) r .
(124) Following [36], and after some computations, we reach the boundary conditions for the kernel equations as follows: It turns out that the observer kernel equation can be transformed into the control kernel equation, therefore obtaining a similar explicit result. For this, defině P n lm (r, ρ) = ρ n−1 r n−1 P n lm (ρ, r), and it can be verified that the equation now governingP n lm (r, ρ) is exactly the equation satisfied by K n lm (r, ρ). ThusP n lm (r, ρ) = K n lm (r, ρ) and we can apply our previous result of Section 5. The observer error dynamics has the same stability properties derived in Section 4 for the closed-loop system under full state control. As in the controller case, only a limited number of modes need to be estimated; namely, those that are not naturally stable by the Lemma 4.1, this being the main difference from the result given in [36].
Finally, the controller-observer augmented system can be proved closed-loop stable as in [36], using the separation principle given the linearity of the system, with desired convergence rate, and without much modification; we skip the details, which requires going up to H 1 stability, as in [36].

Implementation and simulation study
In this section, the simulation experiment on a three-dimensional unity ball (n = 3, R = 1) is taken as an example to illustrate the effectiveness of the proposed control, and some implementation remarks.
The system with the output feedback control law is simulated on 0 ≤ t ≤ 2 s with the following parameters: ε = 1, λ(r) = 10r 4 + 50r 2 + 50, c = 3. We consider that the system initially has the random quantity u 0 ∈ [0, 10], and the observer's initial condition is set as the actual state plus an error of normal distribution with zero mean and σ 2 = 0.5. Fig. 1 shows the plots of the polynomial approximation of kernels K 3 lm , which is obtained by first expanding λ(r)+c ε by using (31), and then finding the coefficients of (32) up to a cutoff in the p-th powers by solving recursively (39)-(41) for each i up to p; in each step one needs to compute the coefficients B i j given by (35) from the previously-found coefficients C i j . The value of K does not depend on m, so we omit this subindex, and l is varied from zero to the value given by the Lemma 3.1. The value of p is chosen as p = 15. Applying the Lemma 3.1, one can obtain l to be 11; however, here, to save space, we only show the first six approximate numerical solutions of the control gains. As shown in Fig. 1, we find that K l becomes increasingly smaller as l increases. In order to avoid a dramatic increase in the complexity of simulation caused by the high dimension, in our simulations, we employ a method also based on spherical harmonic expansions, which greatly reduces the error. Thus, we only calculate the harmonics u m l that only need discretization in the radial direction, and then we sum up a finite number S of harmonics to recover u. When S > 0 is a large enough integer, the error caused by the use of a finite number of harmonics is much smaller than the angular discretization error. Thus, the simulation is carried out using the formula which is a truncated variant of (6), where the spherical harmonics are defined as with P 3 lm the associated Legendre polynomial defined as  Fig. 2(d) and Fig. 3(e), respectively. Note that in these figures, the ranges of color bars are different and thus avoid too uniform colors in Fig. 3. When the open-loop and closed-loop evolution is compared directly, the validity of the proposed method is illustrated more intuitively. Fig. 3(f) shows the L 2 norm of the observation error, from which it can be found that the system begins to converge to its zero equilibrium after the observation error has already settled to zero as well. The evolutions at different layers, namely r = 0.002, r = 0.3, r = 0.5, and r = 0.8, are shown in Fig. 4 (a)-(c), and the observer errors are presented in Fig. 4 (b), (d). For clarity, only the first 0.4 s of the response are shown here. Fig. 5 shows the control effort at the boundary. It can be seen that the system driven by the proposed boundary control eventually converges after a short-term fluctuation.

Conclusion
We have shown a design to stabilize a radially varying reaction-diffusion equation on an n ball, by using an output-feedback boundary control law (with boundary measurements as well) designed through a backstepping method. The radially varying case proves to be a challenge, as the kernel equations become singular in radius; when applying the backstepping method, the same type of singularity appears in the kernel equations, and successive approximations become difficult to use. Using a power series approach, a solution is found, thus providing a numerical method that can be readily applied to both control and observer boundary design. In addition, the required conditions for the radially varying coefficients (analyticity and evenness) are revealed. This result can be extended in several ways. If one has Neumann boundary conditions at the controlled boundary (which implies then that one is measuring the state at the boundary instead of its normal derivative), or even Robin boundary conditions, the method can be extended straightforwardly, since the transformation itself does not change and, therefore, the backstepping kernels remain the same. Only the particular control/observer gains, deduced from the backstepping kernels, would change; as well as a small modification on Lemma 4.1 to account for the change in the boundary conditions.
In practice, this result can be of interest for the deployment of multi-agent systems, following the spirit of [29]; thus, the radial domain mirrors a radial topology of interconnected agents that follow the reaction-diffusion dynamics to converge to equilibria that represent different deployment profiles. Since the plant can be chosen as desired (thereby setting the behavior of the agents), the use of analytic reaction coefficients is not actually a restriction, but opens the door to richer families of deployment profiles compared to the constant coefficient case of [29].
However, the theoretical side of the result needs to be further investigated; an avenue of research that can be explored is the relaxation of the analyticity hypothesis by using reaction coefficients belonging to the Gevrey family; the kernels can then be analyzed to verify if they are still analytic, or rather Gevrey-type kernels, or simply do not converge. Also, the rate of convergence of the obtained power series is of interest and shall be explored. We have experimentally observed that the rate of convergence of the series representation of λ(r) has a considerable influence on it. In addition, one could also explore how fast the series converges in the case of constant λ, since an explicit solution is known from [36]. In particular, the worst case in a domain with radius R would be given by the convergence rate of the Maclaurin series of I 1 λ ε R where I 1 is a modified Bessel function of order 1. Since this function behaves quite closely to an exponential if its argument is large (which would be the case with slowest convergence), the number of required terms would be given by the remainder of the power series of an exponential. In that case, it is easy to see that the size of the term λ ε R would define the required truncation level. If we extrapolate this behavior, then, beyond the convergence rate of the series representation of λ(r), we can say that higher values of R and max r∈[0,R] |λ(r)| and lower values of ε would result in slower-converging series; coincidentally, these are exactly the same factors that would result in a more unstable open-loop plant.