extfieldsofstudy:
Computer Science
Mathematics

text:
Multi-Adaptive Galerkin Methods for ODEs I

We present multi-adaptive versions of the standard continuous and discontinuous Galerkin methods for ODEs. Taking adaptivity one step further, we allow for individual time-steps, order and quadrature, so that in particular each individual component has its own time-step sequence. This paper contains a description of the methods, an analysis of their basic properties, and a posteriori error analysis. In the accompanying paper [A. Logg, SIAM J. Sci. Comput., 27 (2003), pp. 741-758], we present adaptive algorithms for time-stepping and global error control based on the results of the current paper.

1. Introduction. In this paper, we present multi-adaptive Galerkin methods for initial value problems for systems of ODEs of the form where u : [0, T ] → R N , f : R N × (0, T ] → R N is a given bounded function that is Lipschitz-continuous in u, u 0 ∈ R N is a given initial condition, and T > 0 is a given final time. We use the term multi-adaptivity to describe methods with individual time-stepping for the different components u i (t) of the solution vector u(t) = (u i (t)), including (i) time-step length, (ii) order, and (iii) quadrature, all chosen adaptively in a computational feedback process. In the companion paper [29], we apply the multi-adaptive methods to a variety of problems to illustrate the potential of multiadaptivity.
The ODE (1.1) models a very large class of problems, covering many areas of applications. Often different solution components have different time-scales and thus ask for individual time-steps. A prime example to be studied in detail below is our own solar system, where the moon orbits around Earth once every month, whereas the period of Pluto is 250 years. In numerical simulations of the solar system, the time-steps needed to track the orbit of the moon accurately are thus much less than those required for Pluto, the difference in time-scales being roughly a factor 3,000.
Surprisingly, individual time-stepping for ODEs has received little attention in the large literature on numerical methods for ODEs; see, e.g., [4,21,22,3,34]. For specific applications, such as the n-body problem, methods with individual time-stepping have been used-see, e.g., [31,1,5] or [25]-but a general methodology has been lacking. Our aim is to fill this gap. For time-dependent PDEs, in particular for conservation laws of the typeu + f (u) x = 0, attempts have been made to construct methods with individual (locally varying in space) time-steps. Flaherty et al. [20] have constructed a method based on the discontinuous Galerkin method combined with local forward Euler time-stepping. A similar approach is taken in [6], where a method based on the original work by Osher and Sanders [33] is presented for conservation laws in one and two space dimensions. Typically the time-steps used are based on local CFL conditions rather than error estimates for the global error and the methods are low order in time (meaning ≤ 2). We believe that our work on multi-adaptive Galerkin methods (including error estimation and arbitrary order methods) presents a general methodology to individual time-stepping, which will result in efficient integrators also for time-dependent PDEs.
The methods presented in this paper fall within the general framework of adaptive Galerkin methods based on piecewise polynomial approximation (finite element methods) for differential equations, including the continuous Galerkin method cG(q) of order 2q, and the discontinuous Galerkin method dG(q) of order 2q + 1; more precisely, we extend the cG(q) and dG(q) methods to their multi-adaptive analogues mcG(q) and mdG(q). Earlier work on adaptive error control for the cG(q) and dG(q) methods include [7,16,24,18,17,19]. The techniques for error analysis used in these references, developed by Johnson and coworkers (see, e.g., [11,12,10,13,14,15], and [8] in particular) naturally carries over to the multi-adaptive methods.
The outline of the paper is as follows: In section 2 we summarize the key features of the multi-adaptive methods, and in section 3 we discuss the benefits of the new methods in comparison to standard ODE codes. We then motivate and present the formulation of the multi-adaptive methods mcG(q) and mdG(q) in section 4. Basic properties of these methods, such as order, energy conservation, and monotonicity, are discussed in section 5. In the major part of this paper, section 6, we derive a posteriori error estimates for the two methods based on duality arguments, including Galerkin errors, numerical errors, and quadrature errors. We also prove an a posteriori error estimate for stability factors computed from approximate dual solutions.
2. Key features. We summarize the key features of our work on the mcG(q) and mdG(q) methods as follows.
2.1. Individual time-steps and order. To discretize (1.1), we introduce for each component, i = 1, . . . , N , a partition of the time-interval (0, T ] into M i subintervals, I ij = (t i,j−1 , t ij ], j = 1, . . . , M i , and we seek an approximate solution U (t) = (U i (t)) such that U i (t) is a polynomial of degree q ij on every local interval I ij . Each individual component U i (t) thus has its own sequence of time-steps, {k ij } Mi j=1 . The entire collection of individual time-intervals {I ij } may be organized into a sequence of time-slabs, collecting the time-intervals between certain synchronised time-levels common to all components, as illustrated in Figure 2.1.

Global error control.
Our goal is to compute an approximation U (T ) of the exact solution u(T ) at final time T within a given tolerance TOL > 0, using a minimal amount of computational work. This goal includes an aspect of reliability (the error should be less than the tolerance) and an aspect of efficiency (minimal computational work). To measure the error we choose a norm, such as the Euclidean norm · on R N , or more generally some other quantity of interest (see [32]).
The mathematical basis of global error control in · for mcG(q) is an error representation of the form solution U (t), ϕ(t) is the solution of an associated linearized dual problem, and (·, ·) is the R N scalar product.
Using the Galerkin orthogonality, the error representation can be converted into an error bound of the form are stability factors for the different components, depending on the dual solution ϕ(t), and where k i (t) = k ij , q i (t) = q ij for t ∈ I ij . The error bound may take different forms depending on how T 0 (R, ϕ) dt is bounded in terms of R and ϕ. By solving the dual problem numerically, the individual stability factors S i (T ) may be determined approximately, and thus the right-hand side of (2.2) may be evaluated. The adaptive algorithm seeks to satisfy the stopping criterion with maximal time-steps k = (k i (t)).

Iterative methods.
Both mcG(q) and mdG(q) give rise to systems of nonlinear algebraic equations, coupling the values of U (t) over each time-slab. Solving these systems with full Newton may be quite heavy, and we have instead successfully used diagonal Newton methods of more explicit nature.
2.4. Implementation of higher-order methods. We have implemented mcG(q) and mdG(q) in C++ for arbitrary q, which in practice means 2q ≤ 50. The implementation, Tanganyika, is described in more detail in [29] and is publicly (GNU GPL) available for Linux/Unix [30].
2.5. Applications. We have applied mcG(q) and mdG(q) to a variety of problems to illustrate their potential; see [29]. (See also [27] and [26].) In these applications, including the Lorenz system, the solar system, and a number of time-dependent PDE problems, we demonstrate the use of individual time-steps, and for each system we solve the dual problem to collect extensive information about the problems stability features, which can be used for global error control.
3. Comparison with standard ODE codes. Standard ODE codes use timesteps which are variable in time but the same for all components, and the time-steps are adaptively chosen by keeping the "local error" below a given local error tolerance set by the user. The global error connects to the local error through an estimate, corresponding to (2.2), of the form where S is a stability factor. Standard codes do not compute S, which means that the connection between the global error and the local error is left to be determined by the clever user, typically by computing with a couple of different tolerances.
Comparing the adaptive error control of standard ODE codes with the error control presented in this paper and the accompanying paper [29], an essential difference is thus the technique to estimate the global error: either by clever trial-and-error or, as we prefer, by solving the dual problem and computing the stability factors. Both approaches carry extra costs and what is best may be debated; see, e.g., [32] for a comparison.
However, expanding the scope to multi-adaptivity with individual stability factors for the different components, trial-and-error becomes very difficult or impossible, and the methods for adaptive time-stepping and error control presented below based on solving the dual problem seem to bring clear advantages in efficiency and reliability.
For a presentation of the traditional approach to error estimation in ODE codes, we refer to [2], where the following rather pessimistic view is presented: Here we just note that a precise error bound is often unknown and not really needed. We take the opposite view: global error control is always needed and often possible to obtain at a reasonable cost. We hope that multi-adaptivity will bring new life to the discussion on efficient and reliable error control for ODEs.

Multi-adaptive Galerkin.
In this section we present the multi-adaptive Galerkin methods, mcG(q) and mdG(q), based on the discretization presented in section 2.1.
4.1. The mcG(q) method. The mcG(q) method for (1.1) reads as follows: and where P q (I) denotes the linear space of polynomials of degree ≤ q on I. The trial functions in V are thus continuous piecewise polynomials, locally of degree q ij , and the test functions in W are discontinuous piecewise polynomials that are locally of degree q ij − 1.
Noting that the test functions are discontinuous, we can rewrite the global problem (4.1) as a number of successive local problems for each component: For i = 1, . . . , N , j = 1, . . . , M i , find U i | Iij ∈ P qij (I ij ) with U i (t i,j−1 ) given, such that We notice the presence of the vector U (t) = (U 1 (t), . . . , U N (t)) in the local problem for U i (t) on I ij . If thus component U i1 (t) couples to component U i2 (t) through f , this means that in order to solve the local problem for component U i1 (t) we need to know the values of component U i2 (t) and vice versa. The solution is thus implicitly defined by (4.3). Notice also that if we define the residual R of the approximate i.e., the residual is orthogonal to the test space on every local interval. We refer to this as the Galerkin orthogonality for the mcG(q) method.
Making an ansatz for every component U i (t) on every local interval I ij in terms of a nodal basis for P qij (I ij ) (see the appendix), we can rewrite (4.3) as m } q m=1 ⊂ P q−1 (0, 1) are corresponding polynomial weight functions, and τ ij maps Here we assume that the solution is expressed in terms of a nodal basis with the end-points included, so that by the continuity requirement ξ ij0 = ξ i,j−1,qi,j−1 .
Finally, evaluating the integral in (4.5) using nodal quadrature, we obtain a fully discrete scheme in the form of an implicit Runge-Kutta method: For i = 1, . . . , N , j = 1, . . . , M i , find {ξ ijm } qij m=0 , with ξ ij0 given by the continuity requirement, such that for certain weights {w n } (see the appendix). 4.2. The mdG(q) method. The mdG(q) method in local form, corresponding to (4.3), reads as follows: For i = 1, . . . , N , j = 1, . . . , M i , find U i | Iij ∈ P qij (I ij ), such that , and the initial condition is specified for i = 1, . . . , N , by U i (0 − ) = u i (0). On a global level, the trial and test spaces are given by In the same way as for the continuous method, we define the residual R of the approximate solution U as R i (U, t) =U i (t) − f i (U (t), t), defined on the inner of every local interval I ij , and we rewrite (4.7) in the form We refer to this as the Galerkin orthogonality for the mdG(q) method. Notice that this is similar to (4.4) if we extend the integral in (4.4) to include the left end-point of the interval I ij . (The derivative of the discontinuous solution is a Dirac delta function at the end-point.) Making an ansatz for the solution in terms of some nodal basis, we get, as for the continuous method, the following explicit version of (4.7) on every local interval: or, applying nodal quadrature, where the weight functions, the nodal points, and the weights are not the same as for the continuous method.

4.3.
The multi-adaptive mcG(q)-mdG(q) method. The discussion above for the two methods extends naturally to using different methods for different components. Some of the components could therefore be solved for using the mcG(q) method, while for others we use the mdG(q) method. We can even change methods between different intervals.
Although the formulation thus includes adaptive orders and methods, as well as adaptive time-steps, our focus will be mainly on adaptive time-steps.

Choosing basis functions and quadrature.
What remains in order to implement the two methods specified by (4.6) and (4.11) is to choose basis functions and quadrature. For simplicity and efficiency reasons, it is desirable to let the nodal points for the nodal basis coincide with the quadrature points. It turns out that for both methods, the mcG(q) and the mdG(q) methods, this is possible to achieve in a natural way. We thus choose the q + 1 Lobatto quadrature points for the mcG(q) method, i.e., the zeros of xP q (x) − P q−1 (x), where P q is the qth-order Legendre polynomial on the interval; for the mdG(q) method, we choose the Radau quadrature points, i.e., the zeros of P q (x) + P q+1 (x) on the interval (with time reversed so as to include the right end-point). See [28] for a detailed discussion on this subject. The resulting discrete schemes are related to the implicit Runge-Kutta methods referred to as Lobatto and Radau methods; see, e.g., [3].
5. Basic properties of the multi-adaptive Galerkin methods. In this section we examine some basic properties of the multi-adaptive methods, including order, energy conservation, and monotonicity.

5.1.
Order. The standard cG(q) and dG(q) methods are of order 2q and 2q + 1, respectively. The corresponding properties hold for the multi-adaptive methods, i.e., mcG(q) is of order 2q and mdG(q) is of order 2q + 1, assuming that the exact solution u is smooth. We examine this in more detail in subsequent papers.

Energy conservation for mcG(q).
The standard cG(q) method is energy-conserving for Hamiltonian systems. We now prove that also the mcG(q) method has this property, with the natural restriction that we should use the same time-steps for every pair of positions and velocities. We consider a Hamiltonian system, together with initial conditions for x andẋ. Hereẍ is the acceleration, which by Newton's second law is balanced by the force F (x) = −∇ x P (x) for some potential field P . With u = x and v =ẋ we rewrite (5.1) as The total energy E(t) is the sum of the kinetic energy K(t) and the potential energy P (x(t)), Multiplying (5.1) withẋ it is easy to see that energy is conserved for the continuous problem, i.e., E(t) = E(0) for all t ∈ [0, T ]. We now prove the corresponding property for the discrete mcG(q) solution of (5.2). Theorem 5.1. The multi-adaptive continuous Galerkin method conserves energy in the following sense: Let (U, V ) be the mcG(q) solution to (5.2) defined by (4.3). Assume that the same time-steps are used for every pair of positions and corresponding velocities. Then at every synchronized time-levelt, such as, e.g., T , we have 2 and P (t) = P (U (t)). Proof. If every pair of positions and velocities have the same time-step sequence, then we may chooseV as a test function in the equations for U , to get Similarly,U may be chosen as a test function in the equations for V to get and thus K(t) + P (t) = K(0) + P (0). Remark 5.1. Energy conservation requires exact integration of the right-hand side f (or at least that t 0 (U ,V ) dt + (P (t) − P (0)) = 0) but can also be obtained in the case of quadrature; see [23].

Monotonicity.
We shall prove that the mdG(q) method is B-stable (see [3]).
Theorem 5.2. Let U and V be the mdG(q) solutions of (1.1) with initial data U (0 − ) and V (0 − ), respectively, defined by (4.7) on the same partition. If the righthand side f is monotone, i.e., then, at every synchronized time-levelt, such as, e.g., T , we have Proof. Choosing the test function as v = W = U − V in (4.7) for U and V , summing over the local intervals, and subtracting the two equations, we have Noting that The proof is completed noting that the same analysis applies with T replaced by any other synchronized time-levelt. Remark 5.2. The proof extends to the fully discrete scheme, using the positivity of the quadrature weights.

6.
A posteriori error analysis. In this section we prove a posteriori error estimates for the multi-adaptive Galerkin methods, including quadrature and discrete solution errors. Following the procedure outlined in the introduction, we first define the dual linearized problem and then derive a representation formula for the error in terms of the dual and the residual.
6.1. The dual problem. The dual problem comes in two different forms: a continuous and a discrete. For the a posteriori error analysis of this section, we will make use of the continuous dual. The discrete dual problem is used to prove a priori error estimates.
To set up the continuous dual problem, we define, for given functions v 1 (t) and v 2 (t), where * denotes the transpose, and we note that The continuous dual problem is then defined as the following system of ODEs: with data ϕ T and right-hand side g. Choosing the data and right-hand side appropriately, we obtain error estimates for different quantities of the computed solution. We shall assume below that the dual solution has q continuous derivatives (ϕ (qij ) i ∈ C(I ij ) locally on interval I ij ) for the continuous method and q + 1 continuous derivatives (ϕ (qij +1) i ∈ C(I ij ) locally on interval I ij ) for the discontinuous method.

Error representation.
The basis for the error analysis is the following error representation, expressing the error of an approximate solution U (t) in terms of the residual R(U, t) via the dual solution ϕ(t). We stress that the result of the theorem is valid for any piecewise polynomial approximation of the solution to the initial value problem (1.1) and thus in particular the mcG(q) and mdG(q) approximations.
Theorem 6.1. Let U be a piecewise polynomial approximation of the exact solution u of (1.1), and let ϕ be the solution to (6.3) with right-hand side g(t) and initial data ϕ T , and define the residual of the approximate solution U as R( Assume also that U is right-continuous at T . Then the error e = U − u satisfies (6.4) Proof. By the definition of the dual problem, we have using (6.2) Integrating by parts, we get Thus, with L ϕT ,g (e) = (e(T ), ϕ T ) + T 0 (e, g) dt, we have which completes the proof. We now apply this theorem to represent the error in various norms. As before, we let · denote the Euclidean norm on R N and define v L 1 ([0,T ],R n ) = T 0 v dt. Corollary 6.2. If ϕ T = e(T )/ e(T ) and g = 0, then Corollary 6.3. If ϕ T = 0 and g(t) = e(t)/ e(t) , then

Galerkin errors.
To obtain expressions for the Galerkin errors, i.e., the errors of the mcG(q) or mdG(q) approximations, assuming exact quadrature and exact solution of the discrete equations, we use two ingredients: the error representation of Theorem 6.1 and the Galerkin orthogonalities, (4.4) and (4.9). We first prove the following interpolation estimate.
, then there is a constant C q , depending only on q, such that where π [q] f (x) is the qth-order Taylor expansion of f around x 0 = (a+b)/2, k = b−a, and C q = 1/(2 q q!).
Proof. Using Taylor's formula with the remainder in integral form, we have Note that since we allow the polynomial degree to change between different components and between different intervals, the interpolation constant will change in the same way. We thus have C qi = C qi (t) = C qij for t ∈ I ij .
We can now prove a posteriori error estimates for the mcG(q) and mdG(q) methods. The estimates come in a number of different versions. We typically use E 2 or E 3 to adaptively determine the time-steps and E 0 or E 1 to evaluate the error. The quantities E 4 and E 5 may be used for qualitative estimates of error growth. We emphasize that all of the estimates derived in Theorems 6.5 and 6.6 below may be of use in an actual implementation, ranging from the very sharp estimate E 0 containing only local quantities to the more robust estimate E 5 containing only global quantities.
Theorem 6.5. The mcG(q) method satisfies the following estimates: where (6.10) with C q as in Lemma 6.4, k i (t) = k ij , r i (t) = r ij , and s for t ∈ I ij , (6.11) and where π k ϕ is any test space approximation of the dual solution ϕ. Expressions such as C q−1 k q R are defined componentwise, i.e., Proof. Using the error representation of Theorem 6.1 and the Galerkin orthogonality (4.4), noting that the jump terms disappear since U is continuous, we have where π k ϕ is any test space approximation of ϕ. By the triangle inequality, we have Choosing π k ϕ i as in Lemma 6.4 on every interval I ij , we have Continuing, we have , and, finally, As an alternative we can use Cauchy's inequality in a different way. Continuing from E 2 , we have where |R(U, ·)| denotes the vector-valued function with components |R| i = |R i | = |R i (U, ·)|. Noting now that s is the L 2 -projection of |ϕ (q) | onto the piecewise constants on the partition, we have The proof of the estimates for the mdG(q) method is obtained similarly. Since in the discontinuous method the test functions are on every interval of one degree higher order than the test functions in the continuous method, we can choose a better interpolant. Thus, in view of Lemma 6.4, we obtain an extra factor k ij in the error estimates.
Theorem 6.6. The mdG(q) method satisfies the following estimates: where (6.14) and we otherwise use the notation of Theorem 6.5. Proof. As in the proof for the continuous method, we use the error representation of Theorem 6.1 and the Galerkin orthogonality (4.9) to get By Lemma 6.4 we obtain Continuing now in the same way as for the continuous method, we have E 2 ≤ E 3 ≤ E 4 and E 2 ≤ E 5 . Remark 6.1. When evaluating the expressions E 0 or E 1 , the interpolant π k ϕ does not have to be chosen as in Lemma 6.4. This is only a convenient way to obtain the interpolation constant. In section 6.6 below we discuss a more convenient choice of interpolant.
Remark 6.2. If we replace 1 kij Iij |R i | dt by max Iij |R i |, we may replace C q by a smaller constant C ′ q . The value of the constant thus depends on the specific way the residual is measured.

Computational errors.
The error estimates of Theorems 6.5 and 6.6 are based on the Galerkin orthogonalities (4.4) and (4.9). If the corresponding discrete equations are not solved exactly, there will be an additional contribution to the total error. Although Theorem 6.1 is still valid, the first steps in Theorems 6.5 and 6.6 are not. Focusing on the continuous method, the first step in the proof of Theorem 6.5 is the subtraction of a test space interpolant. This is possible, since by the Galerkin orthogonality we have for all test space interpolants π k ϕ. If the residual is no longer orthogonal to the test space, we add and subtract this term to get to the point where the implications of Theorem 6.5 are valid for one of the terms. Assuming now that ϕ varies slowly on each subinterval, we estimate the remaining extra term as follows: (6.16) whereφ is a piecewise constant approximation of ϕ (using, say, the mean values on the local intervals), is a stability factor, and we define the discrete or computational residual as More precise estimates may be used if needed. For the mdG(q) method, the situation is similar with the computational residual now defined as Thus, to estimate the computational error, we evaluate the computational residuals and multiply with the computed stability factors. 6.5. Quadrature errors. We now extend our analysis to take into account also quadrature errors. We denote integrals evaluated by quadrature with˜ . Starting from the error representation as before, we have for the mcG(q) method (6.20) if the quadrature is exact forU v when v is a test function. The first term of this expression was estimated in Theorem 6.5 and the second term is the computational error discussed previously (where˜ denotes that in a real implementation, (6.18) is evaluated using quadrature). The third term is the quadrature error, which may be nonzero even if f is linear, if the time-steps are different for different components. To estimate the quadrature error, notice that (6.21) are the same stability factors as in the estimate for the computational error and is the quadrature residual. A similar estimate holds for the mdG(q) method. We now make a few comments on how to estimate the quadrature residual. The Lobatto quadrature of the mcG(q) method is exact for polynomials of degree less than or equal to 2q − 1, and we have an order 2q estimate for˜ − in terms of f (2q) , and so we make the assumption R Q ij ∝ k 2qij ij . If, instead of using the standard quadrature rule over the interval with quadrature residual R Q0 ij , we divide the interval into 2 m parts and use the quadrature on every interval, summing up the result, we will get a different quadrature residual, namely where we have dropped the ij subindices. Thus, since Thus, by computing the integrals at two or more dyadic levels, we may estimate quadrature residuals and thus the quadrature error.
For the mdG(q) method the only difference is that the basic quadrature rule is one order better, i.e., instead of 2q we have 2q + 1, so that 6.6. Evaluating E G . We now present an approach to estimating the quantity T 0 (R(U, ·), ϕ − π k ϕ) dt by direct evaluation, with ϕ a computed dual solution and π k ϕ a suitably chosen interpolant. In this way we avoid introducing interpolation constants and computing derivatives of the dual. Note, however, that although we do not explicitly compute any derivatives of the dual, the regularity assumed in section 6.1 for the dual solution is still implicitly required for the computed quantities to make sense. Starting now with for the continuous method, we realize that the best possible choice of interpolant, if we want to prevent cancellation, is to choose π k ϕ such that R i (U, ·)(ϕ i − π k ϕ i ) ≥ 0 (or ≤ 0) on every local interval I ij . With such a choice of interpolant, we would have (6.27) The following lemmas give us an idea of how to choose the interpolant. Lemma 6.7. If, for i = 1, . . . , N , , t) and f i is linear or, alternatively, f = f (U (t), t) is linear and all components have the same timesteps and order, then every component R i (U, ·) of the mcG(q) residual is a Legendre polynomial of order q ij on I ij , for j = 1, . . . , M i .
Proof. On every interval I ij the residual component R i (U, ·) is orthogonal to P qij −1 (I ij ). Since the conditions assumed in the statement of the lemma guarantee that the residual is a polynomial of degree q ij on every interval I ij , it is clear that on every such interval it is the q ij th-order Legendre polynomial (or a multiple thereof).
Even if the rather strict conditions of this lemma do not hold, we can say something similar. The following lemma restates this property in terms of approximations of the residual. Lemma 6.8. LetR be the local L 2 -projection of the mcG(q) residual R onto the trial space, i.e.,R i (U, ·)| Iij is the L 2 (I ij )-projection onto P qij (I ij ) of R i (U, ·)| Iij , j = 1, . . . , M i , i = 1, . . . , N . Then everyR i (U, ·)| Iij is a Legendre polynomial of degree q ij . Proof.
for all v ∈ P qij −1 (I ij ), so thatR i (U, ·) is the q ij th-order Legendre polynomial on I ij .
Proof. We can write the qth-order Legendre polynomial on [−1, 1] as P q (x) = 1 q!2 q D q ((x 2 − 1) q ). Thus, integrating by parts, we have is zero at −1 and 1 for l < q. Assume now that f is a polynomial of degree q on [−1, 1] with the property (6.28). Since are linearly independent on [−1, 1] and orthogonal to the Radau polynomial Q q , {Q q (x), (x + 1), (x + 1) 2 , . . . , (x + 1) q } form a basis for P q ([−1, 1]). If then f is orthogonal to the subspace spanned by {(x + 1) p } q p=1 , we must have f = cQ q for some constant c, and the proof is complete.
, t) and f i is linear or, alternatively, f = f (U (t), t) is linear and all components have the same timesteps and order, then every component R i (U, ·) of the mdG(q) residual is a Radau polynomial of order q ij on I ij for j = 1, . . . , M i .
Proof. Note first that by assumption the residual R i (U, ·) is a polynomial of degree q ij on I ij . By the Galerkin orthogonality, we have which holds especially for v(t) = (t − t i,j−1 ) p with p = 1, . . . , q, for which the jump terms disappear. Rescaling to [−1, 1], it follows from Lemma 6.9 that the residual R i (U, ·) must be a Radau polynomial on I ij . Also for the discontinuous method there is a reformulation in terms of approximations of the residual. Lemma 6.11. LetR be the local L 2 -projection of the mdG(q) residual R onto the trial space, i.e.,R i (U, ·)| Iij is the L 2 (I ij )-projection onto P qij (I ij ) of R i (U, ·)| Iij , j = 1, . . . , M i , i = 1, . . . , N . Then everyR i (U, ·)| Iij is a Radau polynomial of degree q ij .
Proof. SinceR i (U, ·) is the L 2 -projection of R i (U, ·) onto P qij (I ij ) on I ij , it follows from the Galerkin orthogonality that for any v(t) = (t − t i,j−1 ) p with 1 ≤ p ≤ q. From Lemma 6.9 it then follows that R i (U, ·) is a Radau polynomial on I ij .
We thus know that the mcG(q) residuals are (in the sense of Lemma 6.8) Legendre polynomials on the local intervals and that the mdG(q) residuals are (in the sense of Lemma 6.11) Radau polynomials. This is illustrated in Figure 6.1.
From this information about the residual, we now choose the interpolant. Assume that the polynomial order of the method on some interval is q for the continuous method. Then the dual should be interpolated by a polynomial of degree q −1, i.e., we have freedom to interpolate at exactly q points. Since a qth-order Legendre polynomial has q zeros on the interval, we may choose to interpolate the dual exactly at those points where the residual is zero. This means that if the dual can be approximated well enough by a polynomial of degree q, the product R i (U, ·)(ϕ i − π k ϕ i ) does not change sign on the interval.
For the discontinuous method, we should interpolate the dual with a polynomial of degree q, i.e., we have freedom to interpolate at exactly q + 1 points. To get rid of the jump terms that are present in the error representation for the discontinuous method, we want to interpolate the dual at the beginning of every interval. This leaves q degrees of freedom. We then choose to interpolate the dual at the q points within the interval where the Radau polynomial is zero.
As a result, we may choose the interpolant in such a way that we have (6.29) with α ij = ±1, for both the mcG(q) method and the mdG(q) method (but the interpolants are different). Notice that the jump terms for the discontinuous method have disappeared.
There is now a simple way to compute the integrals Iij R i (U, ·)(ϕ i − π k ϕ i ) dt. Since the integrands are, in principle, products of two polynomials for which we know the positions of the zeros, the product is a polynomial with known properties. There are then constants C q (which can be computed numerically), depending on the order and the method, such that (6.30) Finally, note that there are "computational" counterparts also for the estimates of type E 3 in Theorems 6.5 and 6.6, namely for the continuous method and similarly for the discontinuous method. 6.7. The total error. The total error is composed of three parts-the Galerkin error, E G , the computational error, E C and the quadrature error, E Q : As an example, choosing estimate E 3 of Theorems 6.5 and 6.6 we have the following (approximate) error estimate for the mcG(q) method: for the mdG(q) method we have (6.34) These estimates containing Galerkin errors, computational errors, and quadrature errors also include numerical round-off errors (included in the computational error).
Modelling errors could also be similarly accounted for since these are closely related to quadrature errors, in that both errors can be seen as arising from integrating the wrong right-hand side.
The true global error may thus be estimated in terms of computable stability factors and residuals. We expect the estimate for the Galerkin error, E G , to be quite sharp, while E C and E Q may be less sharp. Even sharper estimates are obtained using estimates E 0 , E 1 , or E 2 of Theorems 6.5 and 6.6. 6.8. An a posteriori error estimate for the dual. We conclude this section by proving a computable a posteriori error estimate for the dual problem. To compute the stability factors used in the error estimates presented above, we solve the dual problem numerically, and we thus face the problem of estimating the error in the stability factors.
To demonstrate how relative errors of stability factors can be estimated using the same technique as above, we compute the relative error for the stability factor S ϕ (T ), defined as for a computed approximation Φ of the dual solution ϕ.
To estimate the relative error of the stability factor, we use the error representation of Theorem 6.1 to represent the L 1 ([0, T ], R N )-error of Φ in terms of the residual of Φ and the dual of the dual, ω. In [28] we prove the following lemma, from which the estimate follows.
Lemma 6.12. Let ϕ be the dual solution with stability factor S ϕ (t), i.e., with data ϕ(t) = 1 specified at time t, and let ω be the dual of the dual. We then have the following estimate: Theorem 6.13. Let Φ be a continuous approximation of the dual solution with residual R Φ , and assume that S ϕ (t)/S ϕ (T ) is bounded by C on [0, T ]. Then the following estimate holds for the relative error of the stability factor S Φ (T ): and for many problems we may take C = 1.
Proof. By Corollary 6.3, we have an expression for the L 1 ([0, T ], R N )-error of the dual, so that (6.38) With C defined as above it now follows by Lemma 6.12 that and the proof is complete. Remark 6.3. We also have to take into account quadrature errors when evaluating (6.35). This can be done in many ways; see, e.g., [9].
Appendix A. Derivation of the methods. This section contains some details left out of the discussion of section 4.
A.1. The mcG(q) method. To rewrite the local problem in a more explicit form, let {s n } q n=0 be a set of nodal points on [0, 1], with s 0 = 0 and s q = 1. A good choice for the cG(q) method is the Lobatto points of [0, 1]. Now, let τ ij be the linear mapping from the interval I ij to (0, 1], defined by n } q n=0 be the {s n } q n=0 Lagrange basis functions for P q ([0, 1]) on [0, 1], i.e., .
We can then express U i on I ij in the form and choosing the λ To simplify the notation, we drop the ij subindices and assume that the time-interval is [0, k], keeping in mind that, although not visible, all other components are present in f . We thus seek to determine the coefficients {ξ n } q n=1 with ξ 0 given, such that for m = 1, . . . , q we have n (τ (t))λ We explicitly compute the inverseĀ [ [q] mn a n0 = −1.
The weight functions may be computed analytically for small q, and for general q they are easy to compute numerically.
A.2. The mdG(q) method. We now make the same ansatz as for the continuous method, where the difference is that we now have q + 1 degrees of freedom on every interval, since we no longer have the continuity requirement for the trial functions. We make the assumption that the nodal points for the nodal basis functions are chosen so that for m = 0, . . . , q ij , where we use ξ − ij0 to denote ξ i,j−1,qi,j−1 , i.e., the value at the right end-point of the previous interval. To simplify the notation, we drop the subindices ij again and rewrite to [0, k]. We thus seek to determine the coefficients {ξ n } q n=0 such that for m = 0, . . . , q we have  mn ). Then, switching back to the full notation, we have (A.20) [q] mn λ [q] n , m = 0, . . . , q.
As for the continuous method, this may be somewhat simplified. Lemma A.2. For the mdG(q) method, we have q n=0ā [q] mn λ [q] n (0) = 1.
Proof. As in the proof for the mcG(q) method, assume that the interval is [0, 1]. Since the value of the expression is independent of f we can take f = 0. We thus want to prove that if f = 0, then the solution U is constant. By the Galerkin orthogonality, we have The mdG(q) method thus reads as follows: For every local interval I ij , find {ξ ijn } qij n=0 , such that for m = 0, . . . , q ij we have for certain weight functions {w [q] n } q n=0 ⊂ P q (0, 1).


Simulation of Piecewise Hybrid Dynamical Systems in Matlab

,

Introduction
A hybrid dynamical system is a system containing on the same time continuous state variables and event variables in interaction. We find hybrid systems in different fields. We cite robotic systems, chemical systems controlled by vans and pumps, biological systems (growth and division) and nonlinear electronics systems.
Because of interaction between continuous and discrete aspects, the behavior of hybrid systems can be seen as extremely complex. However, this behavior becomes relatively simple for piece-wise affine hybrid dynamical systems that can, in contrast, generate bifurcation and chaos. There are many examples such as power electronics DC-DC converters.
The common power electronics DC-DC converters are the buck converter and the boost converter. They are switching systems with time variant structure [9].
DC-DC converters are widely used in industrial, commercial, residential and aerospace environments. These circuits are typically controlled by PWM (Piece Wise Modulation) or other similar techniques to regulate the tension and the current given to the charges. The controller decides to pass from one configuration to another by considering that transitions occur cyclically or in discrete time. In order to make the analysis possible, most of mathematical treatments use some techniques that are based on averaging or discretization. Averaging can mean to wrong conclusions on operation of a system. Discrete models do not give any information on the state of the system between the sampled instants. In addition, they are difficult to obtain. In fact, in most cases, a pure analytic study is not possible. Another possible approach to analyze these converters can be done via some models of hybrid dynamical systems. DC-DC converters are particularly good candidates for this type of analysis because of their natural hybrid structure. The nature of commutations of these systems makes them strongly nonlinear. They present specific complex phenomena such as fractals structures of bifurcation and chaos.
The study of nonlinear dynamics of DC-DC converters started in 1984 by works of Brockett and Wood [4]. Since then, chaos and nonlinear dynamics in power electronics circuits have attracted different research groups around the world. Different nonlinear phenomena have been observed such as routes to chaos following the period doubling cascade [16], [5], [19], [20] and [23] or quasi-periodic phenomena [6], [7] and [8], besides border collision bifurcations [23] and [2].
Switched circuits behavior is mostly simulated by pure numerical methods where precision step is increased when the system is near a switching condition. Those numerical tools are widely used mainly because of their ease-of-use and their ability to simulate a wide range of circuits including nonlinear, time-variant, and non-autonomous systems.
Even if those simulators can reach the desired relative precision for a continuous trajectory, they can miss a switching condition and then diverge drastically from the trajectory as in figure 1. This could be annoying when one is interested by border collision bifurcations, or when local behavior is needed with a good accuracy. In those applications, an alternative is to write down analytical, or semi-analytical, trajectories and switching conditions to obtain a recurrence which is very accurate and fast to run. Building and adapting such ad'hoc simulators represent a lot of efforts and a risk of mistakes.
Generic and accurate simulators can be proposed if we are restricted to a certain class of systems. A simulation tool with no loss of events is proposed in [14] and [15] for PWAHSs defined on polytopes (finite regions that are bounded by hyperplanes). This class of PWA differential systems has been widely studied as a standard technique to approximate a range of nonlinear systems.
But closed polytopic partition of the state space does not allow simulation of most switching circuits where switching frontiers are mostly single affine constraints or time-dependent periodical events.
We focus on planar PWAHSs with such simple switching conditions which can model a family of switched planar circuits: bang-bang regulators, the Boost converter, the Charge-Pump Phase Locked Loop (CP-PLL), . . . This class of systems has analytical trajectories that help to build fast algorithms with no loss of events. We propose a semi-analytical solver for hybrid systems which provides: • A pure numerical method when the system is nonlinear or non-planar; • A pure analytic method when all continuous parts of the system and switching conditions can be solved symbolically. This can be the case for the boost converter [3], [21], the second order charge-pump phase locked loop [17], [22].
• A mixed method using analytical trajectories and numerical computation of the switching instant when those solutions are transcendent. This has been used for the third order CP-PLL [17]. It can also be the case for the Buck converter [10], [21], . . .
This chapter is organized as follows. Section 2 contains our main results. We describe the problem to be deal with, we introduce a general algorithm to solve planar HSs, we present the algorithm that detects events' occurrence and devote a subsection to our approach efficiency. Section 3 Illustrates the current-mode controlled Boost converter example. Finally, a conclusion is stated in Section 4.

A HS (X, E, t): general definition
A general definition of HSs is presented here. This type of dynamical systems is characterized by the coexistence of two kinds of state vectors: continuous state vector X(t) of real values, and discrete state vector E(t) belonging to a countable discrete set M. Definition 1. A continuous-time, autonomous HS is a system of the form:

HSs class of interest
We consider a two dimensional PWAHS (X(t) ∈ R 2 ). F has then the affine piece-wise form, where A E(t) ∈ R 2× 2 and B E(t) ∈ R 2 are two matrices that depend on the discrete state E(t). Hence, a two dimensional PWAHS is a HS that take the form: We consider two kinds of events: state dependent events and periodic events.
The state dependent event transition S E i E j is defined by an affine state border of the form N ij .X < l ij . In this case an event can occur when the continuous state reaches the border of Note that the set S E i E j is not polytopic in the sense that the domain is not the interior of a closed bounded polytope.

Remark -We consider, without loss of generality, the case where a transition occurs at time dS E i E j if and only if the state X(dS E i E j ) reaches a border of the set S
T is the period and ϕ is the phase of such periodic events.

Event-driven simulation of PWAHSs
The simulation will compute the hybrid state from event to event. Knowing the states X(t k ) and . Then the following algorithm runs the simulation determining the date at the next event as the smallest: ; end Algorithm 1: Algorithm computing the hybrid state at t k+1 .

Event detection occurrence: description and algorithm
We consider the affine Cauchy problem in R 2 : where X 0 is the initial value. We compute the smallest strictly positive time t * i so that the trajectory of X(t) intersects the fixed border B i arriving from the part of the plan where N i .X < l i . The function f i (t) = N i .X(t) − l i defines the guard condition for a border B i . Thus, the problem can be formulated as follows: If f i does not have any strictly positive root or the last condition is not satisfied, t * i is given the infinite value.

Analytical trajectories
Definition 2. For any square matrix A of order n and t in R, the exponential matrix e tA is defined by where I is the identity matrix.
It is well known that the analytical trajectory X(t) of the initial value matrix differential equation (3) is given in terms of the exponential matrix and the variation of constants formula by the general integral form: When A is invertible, the above expression becomes linear: The analytical expression of the exponential matrix e At takes two forms depending on whether the eigenvalues p 1 and p 2 of the matrix A are equal or not: where A = a 11 a 12 a 21 a 22 , A • = a 22 −a 12 −a 21 a 11 and I = 1 0 0 1 .
Using these expressions, we can determine the function f (t) of the problem (4) as follows: where a i are real scalars.
Depending on the eigenvalues p 1 and p 2 , there are five cases that determine the values of the coefficients a i as shown in Table 1. Remark -Coefficients a i are real scalars that depend on Table 1. Expressions of f (t) depending on the eigenvalues p 1 and p 2 .
the eigenvalues p 1 and p 2 , the initial point X k and the border parameters are N i and l i .
In some cases, (p 1 = p 2 = 0, gray cell in Table 1) roots of f (t) can be found analytically and the problem is solved with machine precision.
In other cases, the solution can not be found with classical functions and then a numeric algorithm should be used. Using classical methods like Newton does not guaranty existence or convergence of the smallest positive root. To meet these conditions, let us use analytical roots of the derivative function f (t) expressed in Table 2. We can then compute analytically Table 2. Expressions of f (t) depending on the eigenvalues p 1 and p 2 the set L of ordered roots of f (t), those roots determines monotone intervals of f (t). The following algorithm is used to return the solution t * when it exists or the value ∞ if not.
Remark -When (p 1 , p 2 ) ∈ C * × C * the set L is infinite: when the real part of p i is positive, the algorithm Break; end end Algorithm 2: Algorithm computing t * when a solution is transcendent.
will end by finding a root. In the other case, the set L should be reduced to its three first elements, to find a crossing point when it exists.

Matlab modelling
Our semi-analytical solver is composed of different main programs that define the studied affine system. First, we create the affine system given with a specifically chosen name. Then, we define the matrices A i and B i . After that, we give the switching borders with the sign of transitions and all necessary elements or we give the period if it is about a periodic event. Finally, we execute the simulation by specifying the initial state and the time of simulation.

Application: Current-mode controlled Boost converter
A current-mode controlled Boost converter in open loop consists of two parts: a converter and a switching controller. The basic circuit is given in figure 3.
This converter is a second-order circuit comprising an inductor, a capacitor, a diode, a switch and a load resistance connected in parallel with the capacitor.
The general circuit operation is driven by the switching controller. It compares the inductor current i L with the reference current I re f and generates the on/off driving signal for the switch S. When S is on, the current builds up in the inductor.
When the inductor current i L reaches a reference value, the switch opens and the inductor current flows through the load and the diode. The switch is again closed at the arrival of the next clock pulse from a free running clock of period T.
The Boost converter controlled in current mode is modeled by an affine piece-wise hybrid system defined by the same sub-systems given in equation as follows: In the case of the Boost converter controlled in current mode, there are two types of events: A state event defined by a fixed border of the set S E 1 E 2 : and another periodic event defined by the dates t = dP E i E j , where dP E i E j belongs to the set: where T is the period of this periodic event. The different simulations are obtained using our planar PWA solver.
Before performing any study of the observed bifurcations in this circuit, a numerical simulation in the parametric plane is needed.
The following program calcule_balayage_mod.m is used to obtain the parametric plane:  The figure 4 of the parametric plane allows to emphasize the parameters values for which there exists at least one attractor (fixed point, cycle of order k, strange attractor). Figure 5 shows a bifurcation diagram (Feigenbaum type) in the plane (I re f , i L ). However, figure 6 shows the bifurcation diagram in the space plane (I, i L , v C ).
To draw these two figures we use programs: calcule_figuier.m and affiche_figuier.m   In these two figures, the voltage V in is fixed to 10V and the current I re f varies in the interval [0.5, 1.6]. We observe a period cascade doubling leading to a chaotic regime, interrupted by a border collision bifurcation at I re f = 1.23A (see figure 7). In this figure, a distinction is given between the attractors of attractive cycle type of the order 1 to 14. Each cycle of order k is associated with one color.
For example, the blue area O1 represents the parameters' values for which there exists an attractive fixed point (fundamental periodic regime). The red area O2 represents the existence of an attractive cycle of order 2. The yellow area O4 represents the existence of an attractive cycle of order 4 and so on until getting the cycles O14 of order k = 14. The black area O+ corresponds to parameters values (I re f , V in ) for which there exist cycles of order k ≥ 15 or other types of attractors. In this last area, a chaotic phenomenon could be observed. This bi-dimensional diagram shows some bifurcation curves. In fact, for the rectangle defined by the interval of parameter V in ∈ [7,15] and the parameter I re f ∈ [0.5, 1.6], we observe an area of blue color (existence of attractive fixed point) followed by an area of red color (existence of cycle of order 2), an area of yellow color (existence of cycle of order 4) and another area of black color (existence of cycle of order k ≥ 15 or another attractor type); this succession of zones corresponds to the existence of period doubling cascade.

Conclusion
In this chapter, we have showed an accurate and fast method to determine events' occurrence for planar piece-wise affine hybrid systems. As a result, we have implemented our algorithm in Matlab toolbox version (free downloadable on http://felguezar.000space.com/). This toolbox has also been completed by analysis tools such as displaying the bifurcation and parametric diagrams. The algorithm takes the advantage of the analytical form that appears in the planar case. Our approach can not be extended to a higher dimension. DC-DC converters like Boost converter are known to be simple switched circuits but very rich in nonlinear dynamics. As application, we have chosen the example of Boost converter controlled in current mode.


A novel reduced power compressive sensing technique for wideband cognitive radio

Wideband spectrum sensing for cognitive radio requires high rate analog to digital (A/D) converters whose power consumption is proportional to the sampling rate. In this article, we propose to use sub-Nyquist non-uniform sampling for spectrum sensing to reduce the power consumption. Since the received signal samples are correlated in the time domain, we estimate the missing samples by using the expectation-maximization (EM) algorithm. It is shown that the combined sub-Nyquist non-uniform sampling and EM algorithm consume much less power than A/D converter at the Nyquist rate making the proposed algorithm a viable low-power solution for spectrum sensing. Moreover, it is shown by simulations that the proposed sub-Nyquist rate non-uniform sampler is accurate enough to detect the edges of the estimated power spectral density.

Introduction
A cognitive radio (CR) system [1,2] improves the spectrum utilization by allowing secondary users (SU) to access unused licensed spectrum. Moreover, it assesses its environment and adapts its parameters (e.g. the frequency and power of transmission) to reduce power consumption while having a reliable communication. One reason for increased power consumption is due to electromagnetic pollution. As more wideband access systems are designed, the underlying electromagnetic noise floor increases which in turn requires licensed users to increase their transmit power to maintain required signal to noise ratios (SNRs). CR adds little if any to the noise floor by allowing SU access to the unused bands. Furthermore, by reducing the power consumption of the SU themselves, we can significantly increase transmitted data volumes without a large increase in power consumption or CO 2 emissions. Along with reducing people's exposure to electromagnetic waves, this makes it an excellent candidate for green technology [3]. One important way to reduce power consumption is through developing new signal processing algorithms with lower power consumption. The vacant subbands used by SU are found using spectrum sensing. The first step in spectrum sensing is power spectral density (PSD) estimation. In non-parametric PSD estimation techniques, the PSD is obtained by observing the finite samples of the signal itself. Discrete Fourier transform (DFT)-based PSD estimation methods are the most widely used one [4]. Welch's method [5,6] is an example of DFT-based PSD estimation that was used for spectrum sensing for CR [7]. Simplified DFT (SDFT) is also proposed for PSD estimation in CR [8]. These methods are the simplest methods of PSD estimation, and thus are proposed for PSD estimation in CR [7].
In parametric PSD estimation methods, the received signal samples are modeled as the output of a linear filter to white noise. The PSD is then obtained by estimating the coefficients of that filter. The most important examples are the Yule-Walker auto-regressive (AR) [9] method and the Burg method [9]. These methods are more complex than the previously mentioned non-parametric methods, and thus have not yet been suggested for CR applications.
In subspace methods, the PSD is estimated based on eigenvalue analysis of the correlation matrix. The most important ones are the MUltiple SIgnal Classification (MUSIC) method [9] and the Eigenvector (EV) method [9]. In these methods, first the autocorrelation matrix of http://jwcn.eurasipjournals.com/content/2012/1/281 the received signal is calculated. Then, the EVs of the autocorrelation matrix are computed. The largest eigenvalues correspond to the signal subspace and are used to estimate the signal spectrum. These methods are even more complex than parametric methods.
The second step in spectrum sensing is to detect the PSD edges [10]. PSD edges are the locations of abrupt transitions between occupied and unoccupied subbands. It is shown in [11] that the edges can be detected by detecting the local extrema of the first derivative of the continuous wavelet transform (CWT) of the PSD with respect to frequency. Alternatively, the edges can be detected by differentiating the PSD with respect to frequency. Since the abrupt change in the PSD happens in the edge locations, the first derivative of the PSD exhibits local extremum at the edge locations. Different methods of edge detection for spectrum sensing have been introduced in [10]. In this article, to simplify the analysis as is done in [11], we initially assume that the received signal PSD in CR has a piecewise rectangular-like shape shown in Figure 1. Later, we use a more realistic shape of Sinc square for PSD.
Irrespective of the above steps, wideband spectrum sensing requires high rate analog to digital (A/D) converters with the associated high power consumption. Compressive sensing is a candidate for reduced power consumption in CR [12].
Compressive sensing uses a sub-Nyquist rate sampler to sense the received signal [13]. It was shown in [14] that despite the sub-Nyquist sampling rate, compressive sensing can recover the signals that are sparse or nearly sparse in one domain with some limits introduced in [15]. In [16], a simple least square reconstruction technique was used to recover the samples. Since the frequency edge vector of the received signal in CR has a nearly sparse nature, compressive sensing can be used to recover the edge locations of the PSD in CR [17]. The recovery algorithm of compressive sensing (l 1 minimization algorithm) is computationally complex, and thus time and energy consuming and cannot be easily implemented for real-time applications such as in CR.
In this article, we propose a new algorithm for compressive sensing to be used in spectrum sensing in CR. We show by simulations and analysis that the proposed sub-Nyquist rate non-uniform sampler allows for accurate detection of the edges of PSD and consumes much less power than the conventional sensing method. The proposed sampler samples only some portions of the received signal and switches off the A/D converter based on a predefined pattern to reduce power consumption. Since the received signal samples are correlated, the missing samples can be estimated. We propose to use expectationmaximization (EM) technique [18,19] to estimate these samples. It is shown that the combined sub-Nyquist sampler and EM algorithm consume much less power than Nyquist rate A/D converter making the proposed algorithm a viable low-power solution for spectrum sensing and thus extending the battery life of the CR.

System model
The compressed sampling procedure can be expressed in matrix format as where x c is the compressed sampling received signal vector of size P × 1, x is the Nyquist rate sampled received signal vector of size N × 1 (P N), and β is the compressive sampling matrix of size P × N. When β = I N (I N is the identity matrix of size N), the Nyquist sampling rate is achieved. If β = D N (D N is the DFT matrix of size N), the frequency domain sampling is obtained. For compressive sensing purposes, we can simply eliminate some rows of the identity matrix to achieve the matrix β of the compressed dimension. This elimination may be done on a random basis as is done in [17]. For example, for a 4 × 4 identity matrix with compression ratio of 0.75, we have in which the third row of identity matrix is deleted.

Spectrum sensing
In Welch's method, the received signal samples are segmented into M segments of length N. The estimated PSD (S) is achieved by linearly averaging the periodograms of all segments as where S m is the mth segment's periodogram. The periodogram of each segment is given by where X m is the DFT sequence of the mth segment defined as and we have The estimation error variance of Welch's method is inversely proportional to the number of segments (M) involved in the averaging process. We have [6] in which σ 2 S represents the estimation variance of Welch's method and σ 2 S m is the one for the mth segment's periodogram.
After PSD estimation, the edges are detected. Different methods of edge detection for wideband spectrum sensing are developed in [10]. A modified version of windowaveraging based edge detection technique proposed in [10] has been applied in this article. The proposed algorithm can be summarized as follows (1) The estimated PSD (in dB scale) is averaged over consecutive non-overlapping windows.
(2) The integer part of the window-averaged PSD is calculated. (3) Based on the low spectrum occupancy assumption in CR, the most frequent integer of windowaveraged PSD (the mode) is assumed to be the most probable noise level. The edges are located in the cross section of the noise level and the PSD curve. (4) The area under the PSD curve (the integral of the PSD) between two consecutive edges (detected in the previous step) gives the energy level of that subband. (5) A subband with energy level above the noise level determined in the previous step is assumed to be an occupied subband. Therefore their edges are rough estimates of the edges of the PSD. (6) The exact edge location can be detected by locating the extremum of the first derivative of the PSD within that specified window.
Steps 4 and 5 in the above algorithm are modifications to the window-averaging based edge detection algorithm presented in [10].

EM algorithm
The EM algorithm is an algorithm implementing maximum likelihood estimation. It can be applied to a set of data when their stochastic model is known, although the parameters of the model might be unknown [19]. The EM is an iterative algorithm which works as follows [19]: (1) Calculate the expected value of the log-likelihood function of the conditional probability distribution of the missing variables given the observed ones. This expected value is considered to be the current estimate of the missing points. (2) Maximize the above-mentioned distribution with respect to the parameters of the estimation such as mean, variance and covariance of the log-likelihood function. The parameters maximizing the distribution are used for the next expectation step. (3) Iterate the above-mentioned steps until convergence. (4) The maximum likelihood estimation of the missing variables is obtained after convergence.

Estimation of missing samples in compressive sensing for CR
In this article, we propose to estimate the missing points of the sub-Nyquist sampled received signal in the time domain by applying the EM algorithm. A non-uniform sub-Nyquist rate sampler is applied to the received signal of M × N matrix (recall that M is the number of segments and N is the number of points in each segment.). http://jwcn.eurasipjournals.com/content/2012/1/281 As is shown in Figure 2, the non-uniform sampler works as follows: (1) The sampler provides all samples in the first segment of a set of M 1 segments. (2) For the next (M 1 − 1) segments, the sampler samples only the first N 1 points of each segment and is turned off for the rest of this segment. (3) For the (M 1 + 1) th segment, we restart the process by repeating steps 1 and 2.
The fully sampled first segment and the N 1 points in each subsequent segment are used along with the EM algorithm to estimate the missing samples.
Equivalently, for the (m, n) th element of the received signal matrix, the A/D converter is switched on when either takes the remainder of a division and Fix (.) gives the quotient of a division.
Since the received signal samples are mixed with a Gaussian noise, it is assumed that each segment has a Gaussian distribution. Since the received signal samples are correlated, it is assumed that all segments construct a bi-variate Gaussian distribution with the adjacent segments.

EM algorithm for bi-variate Gaussian distribution
Let x 1 be the first segment (the known segment) and x m be the mth segment of the received signal samples. The first parts of each segment are known and the other parts are unknown and to be estimated by EM algorithm. The estimated missing points converge to their actual value after some iterations in EM algorithm [19].
Since, the received signal samples over consecutive windows are assumed to construct a bi-variate Gaussian distribution, the conditional distribution of missing variables given the observed ones has normal distribution with mean [19] (representing the expected value of missing points of x m+1 in expectation step of EM) given by: where μ m is the mean of the mth segment at each iteration given by σ 2 m,(m+1) is the covariance of the mth segment and its consecutive segment given by The variance of the conditional distribution of missing variables in the mth segment given the observed ones in the (m − 1) th segment (σ 2 m,m.m−1 ) is given by [19] in which σ 2 m,m is given by The above-mentioned procedure to estimate the means, variances and covariances of different segments is the maximization step of EM algorithm.
After a few iterations, the unknown elements of each segment are estimated. Each segment is built based on its preceding segment. The preceding segment is either fully known or partially estimated by the EM algorithm. Therefore, it is not required to wait for all segments to be received. This makes the proposed algorithm suitable for real time applications as it does not require a large buffer or long processing time. The required processing time depends on the processor speed.
After estimating the missing samples, Welch's method is applied to estimate the PSD and the edges.

Analysis of applying EM to spectrum sensing
The frequency bin values of each segment using DFT transformation are given by From Equations (13) and (8), we have (∀k) Since N−1 n=0 e −2πjkn/N = 0 and by applying Equation (13) into (14), Equation (14) can be re-written as In the EM algorithm, the expected value of the variable is chosen as its final value (E (X m+1 (k)) = X m+1 (k)). Thus we have Based on Welch's method, the PSD is obtained by The estimation variance of the EM algorithm-based Welch's method equals that of Welch's method [6]. We have in which σ 2 S represents the estimation variance of Welch's method and σ 2 S m is the one for the mth segment's periodogram.
By substituting Equation (16) into (17), we have (∀k) The term |X 1 (k)| 2 is equivalent to its expected value E |X 1 (k)| 2 in the EM algorithm. Thus Equation (19) can be re-written as (∀k) It can be seen from Equation (20) that the PSD obtained by EM algorithm has a scaling factor (shown in brackets) that is frequency-independent and therefore the PSD shape obtained after the application of the EM algorithm is not distorted compared to the one obtained by sampling above the Nyquist rate.

Power consumption comparison
The idea behind intermittently switching off the A/D converter is to reduce the power consumption and increase the battery life. To compare the power consumption of compressive and non-compressive methods, we calculate the power consumption of the EM algorithm that is additional to reduced power of the compressive sensing algorithm.
The A/D power consumption is linearly proportional to the sampling rate and thus to the number of received samples in each time frame [20]. Moreover for zero-crossing based ADC, the static power consumption is zero [21]. We have, where P c , P nc and P EM are the power consumption of compressive sensing method, non-compressive sensing method and EM algorithm, respectively. C r is the compression ratio and is given by From Equations (8)- (11), the whole procedure to update the missing points based on the EM algorithm requires approximately of 6N × M additions and 4N × M multiplications, a total of 10N × M operations (either http://jwcn.eurasipjournals.com/content/2012/1/281 addition or multiplications) for each iteration. Therefore, power consumption of EM algorithm is given by where P i is the power consumption per instruction, O is the number of operations, and I is the number of iterations. The power efficiency (η) of the system is defined by The system is efficient when η > 1.

Reduced power A/D converter system design example
Power consumption of two cases of compressive and noncompressive sensing is subject to proper design of the system, especially proper selection of ADC and processor.
The proposed system is efficient only if P c is much lower than P nc (η > 1) with the proper design of the system. An example is given in this subsection that is one possible design of the system.
For a system with a compression ratio of C r = 0.1 and power efficiency of (η > 1), we have: Equivalently, we have This will give us the boundary limit on the combination of ADC and processor that a designer can choose.
For example, if we select a 1.2 V 250 mW 14 b 100 MS/s digitally calibrated pipeline ADC in 90 nm CMOS which consumes P nc = 250mW [22], then from (27), we require Table 1 Complexity comparison table of non-compressive  and em-based compressive sensing methods   Method Complexity order For the given example a processor with an efficiency of x = 0.4GOPS/mW in order to achieve a power savings of a factor of 4 (η = 4). From [23], it is shown that processors can have efficiencies up to 17.3 GOPS/mW, therefore the above power savings is achievable with readily available processors.
To verify the applicability of the proposed method, simulation results for the above-mentioned case are presented in the upcoming section.

Complexity comparison of compressive and non-compressive sensing methods
Complexity of non-compressive and EM-based compressive sensing methods are compared in Table 1. The complexity order of Welch's method is given by O(N 3 ) [6] and complexity order of EM algorithm is given by 10INM in which I is the number of iterations. M is the number of received signal segments and N is the number of samples in each segment. For the given example, it can be seen that the complexity of EM-based compressive sensing is of the same order as of the non-compressive sensing technique.

Simulation results
For a wideband signal with 29% of spectrum occupancy, the simulation results of obtaining the PSD using both compressive and non-compressive sensing methods are shown in Figure 3 where the actual frequency shaper filter is shown at the top subplot, the PSD obtained by both compressive and non-compressive sensing methods are shown at the bottom subplot. The SNR of the strongest subbands in these simulations is SNR = 2 dB and the other subband is 3 dB weaker. SNR is calculated individually over various subbands and is defined as the ratio of the signal power in each subband to the noise power in that subband. The simulation results are derived from 100 non-overlapping frames of 1024 samples each using Welch's method. For the compressive sensing method, the whole first segment along with the first 128 points of the other segments are sampled. For M 1 = M = 100 segments and N 1 = 128, the compression ratio is 13%. The remaining unknown data are estimated using the EM algorithm for 15 iterations. It can be seen from Figure 3 that although the PSD for the compressed data has lower values than the non-compressive sensing one, the edges of the PSD are maintained. The power consumption calculations for this case are given in the previous section.  For the signal whose PSD is shown in Figure 3, the simulation results of the root mean square error (RMSE) of the edge frequency bins versus the number of iterations of the EM algorithm (used in the new compressive sensing method) is shown in Figure 4. Based on Figure 4, we use 15 iterations as a reasonable number after which the EM algorithm can be considered to have converged.
For the case considered in Figure 3, the edge detection rate versus SNR curves of both compressive and non-compressive sensing techniques along with the periodogram-based PSD estimation technique, are shown in Figure 5. In periodogram-based PSD estimation technique (whose detection rate is shown in the bottom curve of Figure 5), only the first segment is involved in the PSD estimation (one known segment). It can be seen that the detection rate of the compressive sensing method with 13% of compression ratio, is comparable to that of the non-compressive method and it is much better than the periodogram-based spectrum sensing method making it a viable alternative.
We now consider a more realistic PSD shaper filter as shown in Figure 6. The edge detection rate versus SNR curves of both compressive and non-compressive sensing techniques along with the periodogram-based PSD estimation technique are shown in Figure 7. The SNR of the strongest subbands in these simulations is  The remaining unknown data is estimated using the EM algorithm for 15 iterations. Again, it can be seen that the detection rate of the compressive sensing method is comparable to that of the non-compressive method and it is much better than the periodogrambased spectrum sensing method making it a viable alternative.
In order to examine the effects of roll-off factor of the PSD shaper filter on the detection probability, raised- cosine filter is used as the PSD shaper filter. For a low-pass raised-cosine filter given by in which k represents frequency bin, k 0 represents the cutoff frequency bin of the filter, and β is the roll-off factor (0 ≤ β ≤ 1). The edge happens at frequency bin (1 + β)k 0 . For β = 0, the filter has rectangular-like shape and for β = 1, it has a cosine function shape. For a low-pass raised cosine shaper filter with cut-off frequency bin k 0 = 100, the edge detection rate versus PSD shaper filter roll-off factor curves of both compressive and non-compressive sensing techniques are shown in Figures 8 and 9 for SNR = 0 and SNR = 5 dB, respectively. The simulation results are derived from 100 nonoverlapping frames of 1,024 samples each using Welch's method. For the compressive sensing method, the whole first segment along with the first 128 points of the other segments are available. The compression ratio is 13% for M 1 = M = 100 segments and N 1 = 128. The remaining unknown data is estimated using the EM algorithm for 15 iterations. Again, it can be seen that the detection rate of the compressive sensing method is comparable to that of the non-compressive method. The detection rate decreases by an increase in the roll-off factor of the PSD shaper filter. To understand the effects of different structures on the performance, the detection rate of the following scenarios are compared.
(1) Case 1: M 1 = 100 and N 1 = 128; i.e., starting from the first segment, one segment is known in each set of M 1 = 100 segments and the first N 1 = 128 points of the other segments are known. Its detection is shown in Figure 5. C r = 13%. The signal whose frequency shaper filter is shown at the top subplot of Figure 3 is used for the simulations. The simulation results are derived from 100 non-overlapping frames of 1,024 samples each using Welch's method. The edge detection rate versus SNR curves of the proposed compressive sensing algorithm for all compressive sensing cases and the one for non-compressive sensing are shown in Figure 10.
For above-mentioned cases, it can be seen in Figure 10 that the detection rate increases by both increasing the number of known segments M 1 (case 2) and the number of first known points in each segment N 1 (case 3). However, increasing N 1 results in higher detection rate compared to increasing M 1 under the conditions of the above-mentioned cases because a larger N 1 results in more accurate estimation of the covariance of the consecutive segments and thus results in a more accurate estimation of the missing points.

Conclusion
In this article, we apply a sub-Nyquist non-uniform sampler for spectrum sensing for use in CR. It is shown by simulations and analysis that the proposed sub-Nyquist rate non-uniform sampler is accurate enough to detect the edges of PSD and consumes much less power than the non-compressive sensing method. The proposed sampler samples only some portions of the received signal and switches off the A/D converter based on a predetermined pattern to reduce power consumption. Since the received signal samples in time domain are correlated, we estimate the missing samples using the EM technique. The analysis of applying EM technique to spectrum sensing shows that the locations of PSD edges are maintained after estimating the missing points of the received signal using EM algorithm. In a design example, it is shown that the combined sub-Nyquist sampling/ EM algorithm consumes much less power than Nyquist-based A/D converter thus making the proposed algorithm a viable low-power solution for spectrum sensing. It is shown that although the estimated PSD using the proposed compressive sensing method results in lower values than the non-compressive sensing one, the edges of the PSD are maintained. Since in CR, it is the location of the PSD edges that is important rather than the exact shape of PSD, the proposed compressive sensing method can be used as a low-power solution for A/D conversion in wideband CR.

Methods
Simulations were run in MATLAB. To generate the received data with a specified spectrum, a white Gaussian noise process is passed through spectrum shaping filter. This constructs the frequency domain samples of the cognitive radio received signal with a specified spectrum shaping filter. The time domain samples can be obtained by getting the inverse DFT of the frequency domain samples. Then noise is added to the constructed signal to introduce a specific signal to noise ratio. The PSD can be estimated using Welch's method for both DFT and SDFT techniques. The edges can be obtained based on the proposed edge detection technique. The detection rate of the proposed techniques can be obtained by comparing the estimated edge locations to the actual ones. The simulations were run several times to get an accurate detection rate using Monte-Carlo simulation technique.

Competing interests
The authors declare that they have no competing interests.