Singularity Excitations and Initial Value Problem in Continuous LTI Systems

A modern challenge in electrical engineering education is to keep the math at a sufficient level, with a goal to find an optimal balance between calculus competence and operative skills needed for real-life technical applications. It is not uncommon that some gaps emerge during this quest, which makes it difficult for undergraduates to entirely understand topics related on prior knowledge. This paper aims to draw attention to several important moments concerning total response of Continuous Linear Time Invariant systems, which are superficially or incorrectly explained in many textbooks, and to offer logical arrangement which can be easily understood and accepted by students. The base of discussion relies on classical calculus background, particularly on Picard’s theorem on existence and uniqueness. This theorem is rarely mentioned in signals and systems textbooks. However, mathematical models of many types of signals don’t satisfy the condition for continuity, which can easily produce difficulties in the learning process. It is shown that some reported disagreements and issues related to initial conditions, can be easily cleared out by using the smooth transition from classical calculus to mathematics used in system theory. It is also shown that the classical method is always the primary tool, even for determining the impulse response, while the impulse response is unnecessary or sufficient to determine the total system response, regardless of whether the convolution integral is used.


I. INTRODUCTION
After an introductory discussion about general properties of signals and systems, typical undergraduate signals and systems course accentuate the importance of continuous Linear Time Invariant (LTI) systems for which the input (excitation) and output (response) are related through the constant-coefficient linear differential equations (DE) [1]- [7]. At that point, undergraduates are taught to understand time-domain analysis of DE, as an elementary tool for examination of LTI systems. Later on, discussion is extended to state-space system representations, typically in the Laplace transform domain. Although time-domain analysis looks obsolete compared to techniques based on Laplace transform, it still stands as a valuable prerequisite for understanding other modeling methods and more advanced principles in many engineering disciplines.
The associate editor coordinating the review of this manuscript and approving it for publication was Mira Naftaly .
Together with system modeling using DE, undergraduates are taught that excitations are usually discontinuous, which can easily produce singularities on the right-hand side (RHS) of DE. In addition, they are taught that excitations can also contain singularities. As practice shows, this is the moment when students easily miss the logical connection with basic undergraduate mathematics and accordingly gain a poor understanding of all other topics that follow [8]- [10].
To make things worse, some courses suffer from oversimplified explanations or even inaccuracies in the same or other linked topics, what can additionally reduce student's confidence in acquired knowledge [8]- [12].
Generally, singularities in RHS of DE are a known problem, treated in many different ways, based on various demanding calculus methods. Most of these methods are generally considered to be beyond the abilities of an average undergraduate [8], [13], [14]. Ironically, some simplified concepts from a highly advanced mathematical disciplinetheory of distributions [8], [15], are considered to be the most appropriate for treatment of the problem. Set of these concepts are usually referred as generalized calculus (GC). One of the best introductory texts about singularity signals and GC can be found in [8].
In order to facilitate acceptance of GC with application to time-domain analysis of signals and systems, two main approaches are used: application of classical method which utilizes undergraduate calculus tools with the slight impact of GC, and the convolution method -which is strongly based on impulse response concept followed by a convolution integral. If both of them are taught carelessly in the same course, or in a group of related courses, they can raise more questions than they answer.
The problem is clearly illustrated by two opposite discussions in two distinguished textbooks, concerning total response of LTI system, and interpretation of total response components [2], [4]. In these textbooks, both methods are discussed and evaluated.
According to [4] the convolution method looks more laborious compared to the classical method, but its advantages outweigh the extra work because the classical method has a couple of ''serious'' drawbacks: • Inability to separate the total response into zero-state and zero-input components which are important in the study of the systems.
• The classical method is restricted to post-initial conditions, whereas pre-initial conditions are usually known.
• The classical method is restricted to a certain class of inputs; it cannot be applied to any input. Thus, the classical method is not recommended for the study of continuous/discrete-time signals and systems. [4].
Contrary to [4], recommendation found in [2] (mis)leads to quite the opposite conclusion. Therefore, determination of the total system response by convolution is not demonstrated in [2]. In contrast, the classical method is advocated but explained only for the case where post-initial and pre-initial conditions are the same.
In our paper it is shown that disagreements found in [2] and [4], and issues related to initial conditions, can be easily cleared out by using the natural transition from typical undergraduate mathematics, to the mathematics used in system theory. It is also shown that the classical method is always the primary tool, even for determining the impulse response, while the impulse response is unnecessary or sufficient to determine the total system response, regardless of whether the convolutional integral is used. Discussion in the paper is performed with application to SISO continual systems, while a conclusion for discrete case and MIMO systems can be easily done by analogy.

II. LTI SYSTEM DESCRIPTION BY DE
Even though state-space modeling is the primary tool in the community of control engineers, modeling with high order DE is widely thought to be more natural in some other engineering communities [16]- [18]. Additionally, state-space formulation merely obscures the problems that are discussed in the paper [8], [10].
DE with constant coefficients of N -th order can be used as a model for a wide range of LTI systems: where represents excitation, or system input, y(t) is the system response, or system output, a k and b k are constant coefficients, and D is a linear operator defined as D n = d n /dt n . It is assumed that a 0, , a N , b 0 are not zero, and that N ≥ M , [2], [4]. RHS of (1) can be expressed as f (t) = Q(D)x(t), P(r) is called characteristic polynomial, and P(r) = 0 is called characteristic equation. Roots r 1 , r 2 ,. . . , r N of the characteristic equation are called characteristic roots, or natural frequencies of the system. Equation (1) describes the system in the definition interval (t 1 , t 2 ). Together with N auxiliary conditions in the form of y (k−1) (t 0 ) = Y k−1 , k = 1, . . . N , it uniquely represents a system in terms of excitation only if: • continuity of the RHS is satisfied in (t 1 , t 2 ), and • auxiliary conditions are defined inside (t 1 , t 2 ), i.e. t 1 < t 0 < t 2 . This pair of conditions are known as Picard's theorem (PT) of existence and uniqueness [6], [19]- [24]. From the physical point of view, behavior of the system for t > t 0 joins the past with t < t 0 continuously at t = t 0 , whatever may have been the state in the past [24]. However, in the context of practical dynamical systems where transients are usual phenomena, continuity condition is frequently unsatisfied.

III. RHS WITH DISCONTINUITIES AND TYPES OF AUXILIARY CONDITIONS
Let we assume that excitation x(t) appears at t 0 = 0, and takes a form: where u(t) is Unit step function, Q 1 (D)δ(t), is a linear combination of delta impulse δ(t) and its derivates, and g(t) is a smooth function in (t 1 , t 2 ). It is assumed in the paper that smooth function possesses derivatives of all orders in its domain.
Since delta impulse δ(t) and its derivates are called singularity signals, or simply singularities, Q 1 (D)δ(t) is called singular part, whereas g(t)u(t) is called nonsingular part of x(t). Both nonsingular and singular part can contribute to the occurrence of singularities in f (t). If f (t) contains singularities, y(t) can be discontinuous in origin. Sufficient physical interpretation of this phenomenon can be found in most of the reference textbooks [1]- [10].
By decomposition (t 1 , t 2 ) in three parts, three characteristic cases are distinguished [1]- [6]: 1) If f (t) doesn't contain singularities in origin, postinitial and pre-initial conditions are the same, ε can be strictly 0, i.e. t 0 = 0, and auxiliary conditions are called only initial conditions. Since PT is satisfied in both (t 1 , 0) and (0, t 2 ), determination of system response can be performed using classical calculus methods, typically the method of undetermined coefficients, D operator method, or method of variation of constants [20]- [23]. Origin is automatically included in result due to continuity of y(t).

2) If f (t) does contain singularities in the origin and
pre-initial conditions are defined, i.e.
, and x(t) = 0 inside (t 1 , 0). Therefore classical calculus methods can be applied for the system response determination only for t < 0.

3) If f (t) does contain singularities in the origin, and
post-initial conditions are defined, i.e.
inside (0, t 2 ). Therefore classical calculus methods can be applied for the system response determination, but only for t > 0. In summary, classical methods based on PT can be applied in complete (t 1 , t 2 ) for case a), (0, t 2 ) for case c), and (t 1 , 0) for case b). For cases, b) and c), classical calculus methods are not directly applicable in respect to full interval (t 1 , t 2 ).

IV. CLASSICAL SOLUTION OF (1), SUITABLE FOR THE CASES WHERE PT IS SATISFIED
When PT is satisfied in observing interval, the solution of (1), can be found as a sum of two distinct components: where y h (t) is called homogeneous or complementary solution, and y p (t) is called a particular solution. The homogeneous solution y h (t) is in the form of a linear combination of characteristic modes of the system [1]- [5]. For example, if all characteristic roots are real and different, then where C i , i = 1, 2, . . . , N are arbitrary constants and e r i ·t is called i-th characteristic mode. Frequently, this part of the solution is called a natural response of the system. The particular solution y p (t) is not unique; it can be any function that satisfies (1). It consists of noncharacteristic mode terms y p2 (t), and can contain characteristic mode terms y p1 (t): Arbitrary constants K i , i = 1, 2, . . . , N are closely linked to C i : the sum C i + K i is uniquely determined by (1), and set of auxiliary conditions. This is important but usually overlooked or deliberately omitted fact in university textbooks. The simplest form of particular solution where K i = 0 is usually used, and then particular solution consists entirely of noncharacteristic mode terms, i.e. y p (t) = y p2 (t). This form of a particular solution is frequently called a forced response, and the same form is assumed in the paper. Particular solution y p (t) = y p2 (t) can be found using several classical methods. For standard technical signals, the method of undetermined coefficients, and D operator method are usually suggested. In contrast, the method of variation of constants is recommended only for uncommon types of excitation. Mainly, D operator method, also known as annihilator method [4], [21]- [23] is extremely efficient for typical engineering excitations, it can be easily proved, and it can be easily translated to Laplace transform method. Furthermore, the author's experience is that D operator method is readily accepted by undergraduates.
In order to derive a particular solution, (1) can be rearranged in a form: where rational function H (D) = Q(D)/P(D), with D as an argument, is sometimes called the system operator [9]. When excitation is in the form of complex exponentials, where C is a constant, it is easy to prove that the particular solution is in the form [1]- [4]: One can notice that y p (t) given in (7) consists entirely of noncharacteristic mode terms; therefore, it can be always distinguished from y h (t). Finally, arbitrary constants C i , i = 1, 2, . . . , N are determined, so the (3) satisfies auxiliary conditions. By choosing appropriate σ 0 , and ω 0 in (7), many useful engineering excitations can be easily processed [4], [9]. Also, if C = C(t) is a polynomial, or g(t) is in the resonance with the system, i.e. r 0 ∈ {r 1 , r 2 , . . . , r N }, some simple rearrangements should be done, and (6) can be used for determination of y p (t) [4], [21]- [23].
It is essential to notice that system operator H () has the same algebraic form as a corresponding transfer function in Fourier or Laplace transform domain. The most properties of H () and particular solution y p (t) can be smoothly and transparently forwarded to Fourier and Laplace transform calculus.

V. FORMS OF A TOTAL SYSTEM RESPONSE IN RESPECT TO AUXILIARY CONDITIONS A. TOTAL SYSTEM RESPONSE, CASE A)
For the case a) total response can be solved separately for t > 0 and for t < 0: y(t) = y 1 (t)u(−t) + y 2 (t)u(t). For t > 0 u(t) =1 and RHS of the (1) is equal to f (t) = Q(D)g(t). Then y 2 (t) = y h (t) + y p (t) with arbitrary constants determined, so the y 2 (t) satisfies auxiliary conditions. As explained earlier, for t < 0, RHS of the (1) can be reduced to x(t) = 0, yielding y p (t < 0) = 0. Hence, for t < 0, y 1 (t) = y h (t) with arbitrary constants determined, so the y 1 (t) satisfies auxiliary conditions.

B. TOTAL SYSTEM RESPONSE, CASE B) 1) CONVOLUTION METHOD, CASE B)
For the cases where PT is not satisfied, classical methods that undergraduates have already learned in introductory mathematical courses are not directly applicable. This problem is usually overcome by GC, which introduces concepts of impulse response and convolution [1]- [6].
In almost all related textbooks, pre-initial conditions are widely explained as system energy stored in system states, and they are considered as separate excitations. Therefore, using the superposition principle, the result can be found by separating the total system response into a sum of zero-input and zero-state components, i.e. y(t) = y ZI (t) + y ZS (t). This terminology is widespread in the literature. Still, there are exceptions, for example, in [2] where the zero-state response is called a forced response, and zero-input response is called a natural response.
The zero-input component y ZI (t) satisfies (1) with RHS equal to zero. Since RHS is continual, PT is satisfied in the complete interval (t 1 , t 2 ), and solution can be determined as a classical solution of homogenous DE, with auxiliary conditions equal to pre-initial conditions. For example, if all characteristic roots are real and different, then where arbitrary coefficients C k are resolved in respect to auxiliary conditions. Zero-state component is a response of a system that is initially at rest, to excitation with possible singularities. It is determined in two phases. In the first phase, with the assumption that the system is causal, the impulse response h(t) of the system is found. In the second phase, the impulse response is used for final zero-state calculation.
Determination of post-initial conditions in (9) can be quite complex for N > 3, and can be entirely avoided utilizing auxiliary impulse response h 1 (t) [1], [4]: The auxiliary impulse response is a causal signal and can be described in a form of h 1 (t) = h A (t)u(t) where h A (t) consists VOLUME 8, 2020 of a linear combination of the system's characteristic modes. Therefore, h(t) can be determined by the use of generalized differentiation:

h(t) = Q(D)h 1 (t) = Q(D)(h A (t)u(t)).
However, for N ≥ M generalized differentiation is unnecessary [1], [4]: In respect to PT, expression (10) has same properties as (9): it is singular for t = 0, and for t > 0 it is regular, homogenous, and has the form P(D)h 1 (t) = 0. But contrary to demanding procedure for determination of post-initial conditions required for (9), post-initial conditions for (10) can be calculated very simply: Derivation of (12) is given in the Appendix.
According to standard undergraduate textbooks, the impulse response h(t) of a LTI system contains a complete inputoutput description of the system. If the impulse response is known, zero-state response to any input x(t) can be found using convolution integral, without continuity restrictions required by PE:

2) CONVOLUTION METHOD, CASE B), EXAMPLE
Let we consider a system described by (1) and the following definitions: with excitation in the form x(t) = u(t) − 2δ(t), and pre-initial conditions y (0 − ) = +14, y(0 − ) = −37/8. Obviously, RHS of DE contain singularities in origin, and pre-initial and post-initial conditions are not equal. The first part is determination of the zero-input response: Undetermined constants are resolved using pre-initial conditions: C 1 = −18/8, C 2 = −19/8. The second part is determination of the zero-state response. In the beginning, determination of auxiliary impulse response is performed, by solving (10), with respect to post-initial conditions defined by (12). Therefore, solution for h 1 (t), t > 0 is given by: Since derivative of the auxiliary impulse response is undetermined constants are: Hence Finally, according to (11) and D operator technique [4], [21]- [23], the impulse response is obtained as: In the end, the zero-state response is calculated by use of convolution: Finally,

3) CLASSICAL METHOD, CASE B)
A disadvantage of the classical solution approach quoted in [4] is the restriction to post-initial conditions, due to violation of PT on (t 1 , t 2 ). If post-initial conditions are determined employing GC, for example by impulse matching method [3], [4], then total system response can be determined separately for (0, t 2 ), with respect to post-initial conditions, and for (t 1 , 0), with respect to pre-initial conditions, since in both subintervals PT is satisfied. Because the procedure for calculation of post-initial conditions with respect to pre-initial conditions (or vice versa) can be quite demanding, it is rarely described in textbooks. For example, it is not described in [2], and is explicitly discouraged in [4]. However, by addition of impulse response concept to the classical method, any total response can be easily determined. Again, the superposition principle is a starting point: y(t) = y ZI (t) + y ZS (t). For zero-input response all conditions are the same as in the previous section and given by (13). Contrary, determination of y ZS (t) should be developed since the classical method can not deal with pre-initial conditions [4].
Let we consider excitation x(t), which is in the form of (2). If Q 1 (D) = 0, then the singular part does not exist, and zero-state response can be found by solving equation: P(D)y ZS (t) = Q(D) (g(t)u(t)) , If RHS of (16) contain singularities, post-initial conditions are not zero and should be calculated using some of the available procedures. Furthermore, according to the impulse matching principle, y ZS (t) does not contain impulses when N ≥ M and can be determined using the simplified procedure. Similar to (10), the auxiliary zero-state response can be calculated by solving equation: There are no singularities in RHS of (17), so it is satisfied that y (k−1) . . N , and the classical method can be used. The result is a causal function where y ZS1A (t) is a smooth function. Finally, the zero-state response is determined as Since y ZS (t) does not contain impulses generalized differentiation is unnecessary: Therefore, despite the fact that RHS of DE contains singularities in origin, if the excitation is nonsingular, there is no need for impulse response calculation or generalized differentiation. In this case, total response determination is more efficient than procedure based on convolution integral.
If excitation contains both singular and nonsingular parts as in (2)  [] Using superposition, we can separately solve a set of two auxiliary equations, both with zero pre-initial conditions: Equation (18) is same as (16), whereas (19) is in the form of the impulse response, so it can be derived as (9), i.e. y ZS3 (t) = h(t), see (9). Finally, the total zero-state response is determined as In this case, generalized differentiation is necessary. It is easy to check that calculation effort needed for the described procedure is similar to calculation of convolution integral.

4) CLASSICAL METHOD, CASE B), EXAMPLE
Let we use pre-initial conditions and same excitation as in previous examples. The zero-input response is the same, and the zero-state response is obtained by solving DE with zero pre-initial conditions y ZS (0 − ) = 0, y ZS (0 − ) = 0: Therefore g(t) = 1, Q 1 (D) = −2.
The auxiliary equation in the form of (17) is The equation doesn't have singularities in RHS, pre-initial and post-initial conditions are the same, PT is valid, and the classical method can be used: Using zero post-initial conditions, undetermined constants are resolved: and The second auxiliary equation is in form of (19). Since h(t) is already calculated in (14): Then, zero state response is completed in a form: Finally, the total response is determined as a sum of zero-state and zero-input responses and has the same form as eq (15).

C. TOTAL SYSTEM RESPONSE, CASE C)
When postnatal conditions are defined, and excitation contains singularities, PT is satisfied for (0, t 2 ), and total response can be evaluated using the classical method for t > 0, i.e. y(t) = y h (t)+y p (t). Zero state response can always be determined independently on auxiliary conditions, using classical solution approach, as well as convolution method.
Since the zero-state response is in the form of (5), it can be used as general form of particular solution as in (3), i.e. y ZS (t) = y p (t), whereas coefficients K i are automatically determined. Then y ZI (t) = y h (t) and arbitrary constants C i from (5) are resolved using post-initial conditions.

VI. SEPARATION OF TOTAL RESPONSE INTO PARTS
If all roots of the characteristic equation are in left-half of the complex plane, and if t → ∞ then y h (t) →0 and y(t) → y p (t). Furthermore, if g(t) is not in the resonance with the system, i.e. r 0 / ∈ {r 1 , r 2 , . .., r N }, the particular solution is amplified excitation, with complex gain H (r 0 ). In some cases, y p (t) has quite a clear physical meaning, called steady-state response, whereas y h (t) is called transient response.
In addition to the total response separation into zero-state and zero-input components, separation into transient and steady-state components is equally important. Therefore, it is of interest that any method used for determination of total response supports both types of separation: One of the quoted disadvantages [4] of the classical method is the inability to separate the total response into zero-state and zero-input components. Strictly speaking, the classical method doesn't offer such a separation in classical procedure when post-initial conditions are defined -only separation into y p (t) and y h (t) is available. But using any of described procedures, the zero-state response can be determined regardless of the type of the auxiliary condition, and the zero-input component can be calculated as Therefore, in order to achieve separation into zero-state and zero-input response, both methods is applicable in all three cases, a), b) and c).
In case b) straight forward outcome is separation into y ZI (t) and y ZS (t) parts. If there is an interest to separate the total response into transient and steady-state part, i.e. particular and homogenous part, and if there is no resonance, separation can be performed using visual identification of y p (t), since y p (t) doesn't contain any characteristic modes, and have the same form as excitation. Then the homogenous solution can be derived as: If there is resonance in the system and separation is not obvious, particular solution y p (t) in assumed form can be determined using the classical method, since it doesn't depend upon auxiliary conditions, and (22) can be applied.

VII. CONCLUSION
It is shown that both methods, with some augmentation as presented in the paper, offer the same outcome, and that both methods possess the same possibilities in determination of total system response. Any kind of separation is possible using both methods, without need for recalculation of pre-initial conditions to post-initial conditions, or vice -versa.
Even when convolution integral is used for total response determination, the classical method is indispensable: it is used for determination of h(t). Furthermore, complicated translation of post-initial conditions from pre-initial conditions is not needed for h(t) determination: simplified procedure can be used. The bottom line is that total system response can be resolved without convolution integral, but only by using classical method augmented with impulse response concept and rarely generalized differentiation.
Regarding efficiency, if we make a comparison between convolution method and classical method in the most similar scenario, case b), we can see that both methods require an equal number of steps, with similar complexity, see table 1. In all other situations, the convolution method requires more calculation effort.
Both of the methods have some distinct advantages: For convolution method, the impulse response can be calculated only once, and response to various excitations can be found using convolution -but only if convolution integral converges, what is not generally the case. On the other hand, using D operator method in the classical approach, steady-state response for the case of standard excitations can be easily found -but the complete procedure should be repeated for each excitation.
The third objection to the classical method quoted in [4] is the restriction to a standard (particular class) of excitations. But the same objection can also be applied to convolution method: for some form of x(t) it is difficult to derive y p (t), but for some other form of x(t) it is challenging to derive convolution integral, or integral can't converge. From the pedagogical point of view, it is always better to use as least as needed of new concepts (in this case GC) in order to explain and apply other new concepts: augmented classical method and initial value issues can be very gradually derived from undergraduate calculus. As a tool, the classical method is more straightforward than the convolution method: it relies only on differentiation, what is algebraically more straightforward operation than integration. Therefore, treatment of auxiliary conditions proposed in the paper, and the classical method with suggested modifications, are complete and easy for application, which makes them an excellent alternative to approaches found in existing textbooks.
The simplicity and elegance of both methods disappear when it comes to linear time-variant systems. So, the logical continuation of this topic is the extension to time-varying systems and a Green function method, which will be the subject of future research.

APPENDIX
Calculation of impulse response h(t) with the help of auxiliary impulse response h 1 (t), based on a set of equations (10) and (12), is easy to learn and efficient for application. Derivation of (12) is somewhere conducted from general impulse matching principle [1], [4], whereas in some textbooks it is treated as well-known prior knowledge, for example in [25]. An alternative derivation of (12), based on Step response and principles suggested in the paper, is presented in this appendix.
Let we consider unit-step response s(t) of a system at initial rest, defined through following DE: P(D)s(t) = u(t), s (k−1) (0 − ) = 0, k = 1, . . . N . (23) Because of the absence of singularity in origin, post-initial and pre-initial conditions are equal: According to (10) h 1 (t) = D s(t). Therefore (24) can be written in the form: h