A fresh look at linear ordinary differential equations with constant

0 downloads 0 Views 170KB Size Report
Jun 19, 2015 - coefficient ordinary differential equations of any order based on the ... Using this, one finds a particular solution when the forcing term is a polynomial, ... The problem with this approach is that one needs to 'digest' the theory of linear ... the characteristic polynomial, and let λ1,...,λn ∈ C be the roots of p(λ) ...
This article was downloaded by: [Roberto Camporesi] On: 19 June 2015, At: 10:08 Publisher: Taylor & Francis Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

International Journal of Mathematical Education in Science and Technology Publication details, including instructions for authors and subscription information: http://www.tandfonline.com/loi/tmes20

A fresh look at linear ordinary differential equations with constant coefficients. Revisiting the impulsive response method using factorization Roberto Camporesi

a

a

Dipartimento di Scienze Matematiche, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy Published online: 19 Jun 2015.

Click for updates To cite this article: Roberto Camporesi (2015): A fresh look at linear ordinary differential equations with constant coefficients. Revisiting the impulsive response method using factorization, International Journal of Mathematical Education in Science and Technology, DOI: 10.1080/0020739X.2015.1057246 To link to this article: http://dx.doi.org/10.1080/0020739X.2015.1057246

PLEASE SCROLL DOWN FOR ARTICLE Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content. This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms &

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

Conditions of access and use can be found at http://www.tandfonline.com/page/termsand-conditions

International Journal of Mathematical Education in Science and Technology, 2015 http://dx.doi.org/10.1080/0020739X.2015.1057246

A fresh look at linear ordinary differential equations with constant coefficients. Revisiting the impulsive response method using factorization Roberto Camporesi∗ Dipartimento di Scienze Matematiche, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

(Received 3 September 2014) We present an approach to the impulsive response method for solving linear constantcoefficient ordinary differential equations of any order based on the factorization of the differential operator. The approach is elementary, we only assume a basic knowledge of calculus and linear algebra. In particular, we avoid the use of distribution theory, as well as of the other more advanced approaches: Laplace transform, linear systems, the general theory of linear equations with variable coefficients and variation of parameters. The approach presented here can be used in a first course on differential equations for science and engineering majors. Keywords: ordinary differential equations; Green functions; instructional exposition

1. Introduction In view of their many applications in various fields of science and engineering, linear constant-coefficient differential equations constitute an important chapter in the theory and application of ordinary differential equations. In introductory courses on differential equations, the treatment of second or higher order non-homogeneous equations is usually limited to illustrating the method of undetermined coefficients. Using this, one finds a particular solution when the forcing term is a polynomial, an exponential, a sine or a cosine, or a product of terms of this kind. It is well known that the impulsive response method gives an explicit formula for a particular solution in the more general case in which the forcing term is an arbitrary continuous function. This method is generally regarded as too difficult to implement in a first course on differential equations. Students become aware of it only later, as an application of the theory of the Laplace transform [1] or of distribution theory.[2] An alternative approach which is sometimes used consists in developing the theory of linear systems first, considering then linear equations of order n as a particular case of this theory. The problem with this approach is that one needs to ‘digest’ the theory of linear systems, with all the issues related to the diagonalization of matrices and the Jordan form ([3], chapter 3). Another approach is by the general theory of linear equations with variable coefficients, with the notion of Wronskian and the method of the variation of constants. This approach can be implemented also in the case of constant coefficients ([4], chapter 2). However, ∗

Email: [email protected]

 C 2015 Taylor & Francis

2

R. Camporesi

in introductory courses, the variation of constants method is often limited to first-order equations. Moreover, this method may be very long to implement in specific calculations, even for rather simple equations. Finally, within this approach, the occurrence of the particular solution as a convolution integral is rather indirect, and appears only at the end of the theory (see, for example, [4], exercise 4, p. 89). The purpose of these notes is to give an elementary presentation of the impulsive response method using only very basic results from linear algebra and calculus in one or many variables. We write the nth order equation in the form

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

y (n) + a1 y (n−1) + a2 y (n−2) + · · · + an−1 y  + an y = f (x),

(1.1)

where we use the notation y(x) in place of x(t) or y(t), y(k) denotes as usual the derivative of order k of y, a1 , a2 , . . . , an are complex constants and the forcing term f is a complex-valued continuous function in an interval I. When f = 0, the equation is called non-homogeneous. When f = 0, we get the associated homogeneous equation y (n) + a1 y (n−1) + a2 y (n−2) + · · · + an−1 y  + an y = 0.

(1.2)

The basic tool of our investigation is the so-called impulsive response. This is the function defined as follows. Let p(λ) = λn + a1 λn − 1 + ··· + an − 1 λ + an (λ ∈ C) be the characteristic polynomial, and let λ1 , . . . , λn ∈ C be the roots of p(λ) (not necessarily distinct, each counted with its multiplicity). We define the impulsive response g = gλ1 ···λn recursively by the following formulas: for n = 1, we set gλ1 (x) = eλ1 x , for n ≥ 2 we set  gλ1 ···λn (x) = e

x

λn x

e−λn t gλ1 ···λn−1 (t) dt

(x ∈ R).

(1.3)

0

It turns out that g solves the homogeneous Equation (1.2) with the initial conditions y(0) = y  (0) = · · · = y (n−2) (0) = 0,

y (n−1) (0) = 1.

Moreover, the impulsive response allows one to solve the non-homogeneous equation with an arbitrary continuous forcing term and with arbitrary initial conditions. Indeed, if g denotes the impulsive response of order n and if 0 ∈ I, we shall see that the general solution of (1.1) in the interval I can be written as y = yp + yh , where the function yp is given by the convolution integral  yp (x) =

x

g(x − t)f (t) dt,

(1.4)

0 (k)

and solves (1.1) with trivial initial conditions at the point x = 0 (i.e., yp (0) = 0 for k = 0, 1, . . . , n − 1), whereas the function yh (x) =

n−1  k=0

ck g (k) (x)

(1.5)

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

International Journal of Mathematical Education in Science and Technology

3

gives the general solution of the associated homogeneous Equation (1.2) as the coefficients ck vary in C. In other words, the n functions g, g , g , . . . , g(n − 1) are linearly independent solutions of this equation and form a basis of the vector space of its solutions. For n ≥ 2, we will not assume, a priori, existence, uniqueness, and extendability of the solutions of a linear initial value problem (homogeneous or not). We shall rather prove these facts directly, in the case of constant coefficients, by obtaining explicit formulas for the solutions. The proof of (1.4) that we give is by induction on n and is based on the factorization of the differential operator acting on y in (1.1) into first-order factors, along with the formula for solving first-order linear equations. It requires, moreover, the interchange of the order of integration in a double integral; that is, the Fubini theorem. The proof is constructive in that it directly produces the particular solution yp as a convolution integral between the impulsive response g and the forcing term f. In particular, if we take f = 0, we get that the unique solution of the homogeneous initial value problem with all vanishing initial data is the zero function. This implies immediately the uniqueness for the initial value problem (homogeneous or not) with arbitrary initial data. We remark that in the usual approach to linear equations with constant coefficients, the proof of uniqueness relies on an estimate for the rate of growth of any solution y of the homogeneous equation, together with its derivatives y , y , . . . , y(n − 1) , in terms of the coefficients a1 , a2 , . . . , an in (1.2). (See, e.g., [4], Theorems 13 and 14, chapter 2.) As regards (1.5), one can give a similar proof by induction on n using factorization. We shall rather prove (1.5) by computing the linear relationship between the coefficients (k) ck and the initial data bk = y (k) (0) = yh (0) for any given n. If x0 is any point of I, we x x can just replace 0 with x0 in (1.4), and g(k) (x) with g(k) (x − x0 ) in (1.5). The function (k)

yp satisfies then yp (x0 ) = 0 for 0 ≤ k ≤ n − 1, and the relationship between ck and (k) bk = y (k) (x0 ) = yh (x0 ) remains the same as before. We will then look for explicit formulas for the impulsive response g. For n = 1, g is an exponential, for n = 2, 3 it is easily computed. For generic n, one can use the recursive formula (1.3) to compute g at any prescribed order. A simple formula is obtained if the roots of p(λ) are all equal or all different. In the general case, one can prove by induction on k that if λ1 , . . . , λk are the distinct roots of p(λ), of multiplicities m1 , . . . , mk , then there exist polynomials G1 , . . . , Gk , of degrees m1 − 1, . . . , mk − 1, such that

g(x) =

k 

Gj (x)eλj x .

(1.6)

j =1

An explicit formula for the polynomials Gj will be given in the Appendix. We also give the formula for the polynomials Gj based on the partial fraction expansion of 1/p(λ). This formula can be proved by induction, but it is most easily proved using the Laplace transform or distribution theory. Using (1.6) in (1.5) gives the well-known formula for the general solution of (1.2) in terms of the complex exponentials eλj x . The case of real coefficients follows easily from this result. Using (1.6) in (1.4) gives a simple proof of the method of undetermined coefficients for the case in which f (x) = P (x)eλ0 x , where P is a complex polynomial and λ0 ∈ C. This paper is an outgrowth of [5], where the method of factorization was illustrated for second-order equations. Here, we treat equations of arbitrary order n. Since the methodology used here is the same as that of [5], a reading of [5] enables the understanding of the current

4

R. Camporesi

paper. We also refer to [6] and the references therein, for the approach by factorization to linear constant-coefficient differential equations. 2. Linear equations of order n Consider the linear constant-coefficient non-homogeneous differential equation of order n y (n) + a1 y (n−1) + a2 y (n−2) + · · · + an−1 y  + an y = f (x),

(2.1)

k

where y (k) = ddxyk , a1 , a2 , . . . , an are complex constants, and the forcing term f is a complexvalued continuous function in the interval I, i.e., f ∈ C0 (I). When f = 0, we obtain the associated homogeneous equation

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

y (n) + a1 y (n−1) + a2 y (n−2) + · · · + an−1 y  + an y = 0.

(2.2)

We will write (2.1) and (2.2) in operator form as Ly = f(x) and Ly = 0, where L is the linear differential operator of order n with constant coefficients defined by Ly = y (n) + a1 y (n−1) + a2 y (n−2) + · · · + an−1 y  + an y for any function y at least n times differentiable. Denoting by we have  L=

d dx



n + a1

d dx

d dx

n−1 + · · · + an−1

the differentiation operator,

d + an . dx

(2.3)

L defines a map C n (R) → C 0 (R) that to each complex-valued function y at least n-times differentiable over R with continuous nth derivative associates the continuous function Ly. The fundamental property of L is its linearity; that is, L(c1 y1 + c2 y2 ) = c1 Ly1 + c2 Ly2 ,

(2.4)

∀c1 , c2 ∈ C, ∀y1 , y2 ∈ C n (R). This formula implies some important facts. First, if y1 and y2 are any two solutions of the homogeneous equation on R, then any linear combination c1 y1 + c2 y2 is also a solution of this equation on R. In other words, the set V = {y ∈ C n (R) : Ly = 0} is a complex vector space. We will prove that dim V = n, and will see how to get a basis of V. Note that if y ∈ V, then y ∈ C ∞ (R), as follows from (2.2) by isolating y(n) and differentiating any number of times. Second, if y1 and y2 are two solutions of (2.1) on a given interval I ⊂I, then their difference y1 − y2 solves (2.2). It follows that if we know a particular solution yp of the non-homogeneous equation on an interval I , then any other solution of (2.1) on I is given by yp + yh , where yh is a solution of the associated homogeneous equation. We will prove that the solutions of (2.1) are defined on the whole of the interval I in which f is continuous (and are of course of class Cn there). The fact that L has constant coefficients allows one to find explicit formulas for the d λx e = λeλx , solutions of (2.1)–(2.2). For the moment, we only observe (by the formula dx

International Journal of Mathematical Education in Science and Technology

5

∀λ ∈ C) that the complex exponential eλx is a solution of (2.2) if and only if λ is a root of the characteristic polynomial p(λ) = λn + a1 λn−1 + · · · + an−1 λ + an .

(2.5)

Let λ1 , λ2 , . . . , λn ∈ C be the roots of p(λ), not necessarily all distinct, each counted with its multiplicity. The polynomial p(λ) factors as p(λ) = (λ − λ1 )(λ − λ2 ) · · · (λ − λn ). The differential operator L can be factored in a similar way as

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

 L=

d − λ1 dx



   d d − λ2 · · · − λn , dx dx

(2.6)

where the order of composition of the factors is unimportant as they all commute. For example, for n = 2 we have Ly = y + a1 y + a2 y, and 

d − λ1 dx



   d d − λ2 y = − λ1 (y  − λ2 y) dx dx = y  − (λ1 + λ2 )y  + λ1 λ2 y,

which coincides with Ly since λ1 + λ2 = −a1 and λ1 λ2 = a2 . Moreover, it is clear that 

d − λ1 dx



d − λ2 dx



 =

d − λ2 dx



 d − λ1 . dx

The idea is now to use (2.6) to reduce the problem to first-order differential equations. The following result gives a particular solution of (2.1). Theorem 2.1: Let f ∈ C0 (I), and suppose that 0 ∈ I. Let λ1 , λ2 , . . . , λn be n complex numbers (not necessarily all distinct), and let L be the differential operator (2.6). Then, the initial value problem 

Ly = f (x) y(0) = y  (0) = · · · = y (n−1) (0) = 0

(2.7)

has a unique solution, defined on the whole of I, and given by the formula 

x

y(x) =

g(x − t)f (t) dt

(x ∈ I ),

(2.8)

0

where g = gλ1 ···λn is the function defined recursively as follows: for n = 1, we set gλ (x) = eλx (λ ∈ C), for n ≥ 2 we set  gλ1 ···λn (x) = 0

x

gλn (x − t) gλ1 ···λn−1 (t) dt

(x ∈ R).

(2.9)

6

R. Camporesi

The function gλ1 ···λn is of class C∞ on the whole of R. In particular, if we take f = 0, we get that the unique solution of the homogeneous problem 

Ly = 0 y(0) = y  (0) = · · · = y (n−1) (0) = 0

(2.10)

is the zero function y = 0. Proof: We proceed by induction on n. First, we prove that gλ1 ···λn ∈ C ∞ (R). Indeed, this holds for n = 1. Suppose it holds for n − 1. Then, it holds also for n, since by (2.9)  gλ1 ···λn (x) = e

x

λn x

e−λn t gλ1 ···λn−1 (t) dt,

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

0

so gλ1 ···λn is the product of two functions of class C∞ on R. (The integral function of a C∞ function is C∞ .) We now prove (2.8). The result holds for n = 1, in view of the formula  y(x) =

x

eλ(x−t) f (t) dt

(2.11)

0

for the unique solution to the first-order initial value problem 

y  − λy = f (x) y(0) = 0.

Assuming the result holds for n − 1, let us prove it for n. Consider then the problem (2.7) with L given by (2.6) 

 d d − λ1 dx − λ2 · · · dx − λn y = f (x) y(0) = y  (0) = · · · = y (n−1) (0) = 0.

Letting h =



d dx

d dx

− λn y, we see that y solves the problem (2.7) if and only if h solves 

 d d − λ1 dx − λ2 · · · dx − λn−1 h = f (x) h(0) = h (0) = · · · = h(n−2) (0) = 0, d dx

and y solves 

y  − λn y = h(x) y(0) = 0.

(The initial conditions for h = y − λn y follow by computing h , h , . . . , h(n − 2) and setting x = 0.) By the inductive hypothesis and by (2.11), we have  h(x) = 0

x

gλ1 ···λn−1 (x − t)f (t) dt,

International Journal of Mathematical Education in Science and Technology

 y(x) =

7

x

eλn (x−t) h(t) dt. 0

Substituting h from the first formula into the second, we obtain y(x) as a repeated integral (∀x ∈ I): 

x

y(x) =

 gλn (x − t)

0

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

t

 gλ1 ···λn−1 (t − s)f (s) ds dt.

(2.12)

0

To fix ideas, suppose x > 0. Then, in the integral with respect to t we have 0 ≤ t ≤ x, whereas in the integral with respect to s we have 0 ≤ s ≤ t. We can then rewrite y(x) as a double integral  y(x) =

gλn (x − t)gλ1 ···λn−1 (t − s) f (s) ds dt, Tx

where Tx is the triangle in the (s, t) plane defined by 0 ≤ s ≤ t ≤ x, with vertices at the points (0, 0), (0, x), (x, x). In (2.12), we first integrate with respect to s and then with respect to t. Since the triangle Tx is convex both horizontally and vertically, and since the integrand function F (s, t) = gλn (x − t) gλ1 ···λn−1 (t − s) f (s) is continuous in Tx , we can interchange the order of integration and integrate with respect to t first. Given s (between 0 and x) the variable t in Tx varies between s and x, see the picture below. t s

t

x s t x

x

0

s

We thus obtain 

x

y(x) = 0

 s

x

 gλn (x − t)gλ1 ···λn−1 (t − s) dt f (s) ds.

8

R. Camporesi

By substituting t with t + s in the integral with respect to t, we obtain 

x

y(x) = 



0

 gλn (x − s − t)gλ1 ···λn−1 (t) dt f (s) ds

0 x

=

x−s

gλ1 ···λn (x − s)f (s) ds,

0

which is (2.8). For x < 0, we can reason in a similar way and we get the same result.



By induction, one can verify that the function gλ1 ···λn solves the following homogeneous initial value problem:

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015



Ly = 0 y(0) = y  (0) = · · · = y (n−2) (0) = 0,

y (n−1) (0) = 1.

(2.13)

Conversely, if y solves (2.13), then y is given by (2.9). In other words, the solution of the initial value problem (2.13) is unique and is calculable for n ≥ 2 by the recursive formula (2.9). This is proved by induction, reasoning as in the proof of Theorem 2.1. For example, take n = 2 and suppose y solves the problem 

 d − λ1 dx − λ2 y = 0 y(0) = 0, y  (0) = 1.

Letting h =



d dx

d dx

− λ2 y, the function h solves 

− λ1 h = 0 h(0) = 1, d dx

whence h = gλ1 . Therefore, y solves 

y  − λ2 y = gλ1 (x) y(0) = 0,

and by (2.11), we get 

x

y(x) =

gλ2 (x − t) gλ1 (t) dt = gλ1 λ2 (x).

0

The function g = gλ1 ···λn is called the impulsive response of the differential operator L. The name comes from the initial conditions in (2.13) for the case n = 2, namely y(0) = 0, y (0) = 1. We will see later how to compute it explicitly using (2.9). Note that as a function of λ1 , . . . , λn , gλ1 ···λn is symmetric in the interchange λi ↔λj , ∀i, j = 1, . . . , n. This follows from the fact that gλ1 ···λn solves the problem (2.13) with L symmetric in λi ↔λj , ∀i, j. It is interesting to verify directly that the function y given by (2.8) solves the nonhomogeneous problem (2.7). Indeed, using repeatedly the formula d dx

 0

x



x

F (x, t) dt = F (x, x) + 0

∂F (x, t) dt, ∂x

International Journal of Mathematical Education in Science and Technology

9

we get from (2.8) 



x

y (x) = g(0)f (x) + 







0





g  (x − t)f (t) dt

0

x

y (x) = g (0)f (x) +

x

g (x − t)f (t) dt = 

g (x − t)f (t) dt =

0

x

g  (x − t)f (t) dt

0

.. .   x x (n−1) (n−2) (x) = g (0)f (x) + g (n−1) (x − t)f (t) dt = g (n−1) (x − t)f (t) dt y 0 0  x  x y (n) (x) = g (n−1) (0)f (x) + g (n) (x − t)f (t) dt = f (x) + g (n) (x − t)f (t) dt,

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

0

0

where we used (2.13). It follows that 

x

(Ly)(x) = f (x) +

(Lg)(x − t)f (t) dt = f (x),

∀x ∈ I,

0

being Lg = 0. The initial conditions y(k) (0) = 0 for 0 ≤ k ≤ n − 1 are immediately verified from these formulas. For the initial value problem with arbitrary (complex) initial data at x = 0, we have the following result. Theorem 2.2: Let f ∈ C0 (I), 0 ∈ I, and let b0 , b1 , . . . , bn − 1 be n complex numbers. Let L be the differential operator in (2.3), (2.6), and let g be the impulsive response of L. Then, the initial value problem 

Ly = f (x) y(0) = b0 , y  (0) = b1 , . . . , y (n−1) (0) = bn−1

(2.14)

has a unique solution, defined on the whole of I, and given for any x ∈ I by  y(x) =

x

g(x − t)f (t) dt + c0 g(x) + c1 g  (x) + · · · + cn−1 g (n−1) (x),

(2.15)

0

where ⎧ c0 = bn−1 + a1 bn−2 + · · · + an−2 b1 + an−1 b0 ⎪ ⎪ ⎪ ⎪ c ⎪ 1 = bn−2 + a1 bn−3 + · · · + an−3 b1 + an−2 b0 ⎪ ⎪ ⎨ .. . ⎪ = b2 + a1 b1 + a2 b0 c ⎪ n−3 ⎪ ⎪ ⎪ c = b1 + a1 b0 ⎪ ⎪ ⎩ n−2 cn−1 = b0 .

(2.16)

In particular (taking f = 0), the solution of the homogeneous problem 

Ly = 0 y(0) = b0 , y  (0) = b1 , . . . , y (n−1) (0) = bn−1

(2.17)

10

R. Camporesi

is unique, of class C∞ on the whole of R, and is given by yh (x) =

n−1 

(x ∈ R).

ck g (k) (x)

(2.18)

k=0

Proof: The uniqueness follows from the fact that if y1 , y2 both solve the problem (2.14), then their difference y˜ = y1 − y2 solves the problem (2.10), whence y˜ = 0 by Theorem 2.1. To show that y(x) given by (2.15) satisfies Ly(x) = f(x), just observe that any derivative g(k) of g satisfies the homogeneous equation, since 

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

Lg (k) = L

d dx



k g=

d dx

k Lg = 0.

By Theorem 2.1 and the linearity of L, it follows that (Ly)(x) = f (x) + L

 n−1 

 ck g

(k)

k=0

(x) = f (x) +

n−1 

ck (Lg (k) )(x) = f (x).

k=0

There remain to prove the relations (2.16) between the coefficients ck and bk = y(k) (0) (0 ≤ k ≤ n − 1). By computing y , y , . . . , y(n − 1) from (2.15) and imposing the initial conditions in (2.14), we obtain the linear system ⎧ ⎪ ⎪ cn−1 = b0 (n) ⎪ ⎪ ⎪ ⎪ cn−2 + cn−1 g (n) (0) = b1 (n+1) ⎪ ⎨ cn−3 + cn−2 g (0) + cn−1 g (0) = b2 (n) (n+1) c + c g (0) + c g (0) + cn−1 g (n+2) (0) = b3 n−4 n−3 n−2 ⎪ ⎪ ⎪ . ⎪ .. ⎪ ⎪ ⎪ ⎩ (n) c0 + c1 g (0) + c2 g (n+1) (0) + · · · + cn−1 g (2n−2) (0) = bn−1 ,

(2.19)

where we used g(k) (0) = 0 for 0 ≤ k ≤ n − 2, g(n − 1) (0) = 1. The derivatives g(k) (0) with k ≥ n can be computed from the differential equation Lg = 0 by isolating g(n) (x), differentiating and letting x = 0. One gets g (n) (0) = −a1 g (n−1) (0) = −a1 g (n+1) (0) = −a1 g (n) (0) − a2 g (n−1) (0) = a12 − a2 g (n+2) (0) = −a1 g (n+1) (0) − a2 g (n) (0) − a3 g (n−1) (0) = −a13 + 2a1 a2 − a3 , and so on. These relations, substituted in the system (2.19), allow one to solve it recursively and we get exactly the relations (2.16).  Remark 1: By imposing the initial conditions at an arbitrary point x0 of the interval I, we get that the unique solution to the problem 

Ly = f (x) y(x0 ) = b0 , y  (x0 ) = b1 , . . . , y (n−1) (x0 ) = bn−1

International Journal of Mathematical Education in Science and Technology

11

is given, ∀x ∈ I, by  y(x) =

x

g(x − t)f (t) dt +

x0

n−1 

ck g (k) (x − x0 ),

(2.20)

k=0

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

where the coefficients ck are given by (2.16). The proof of this is analogous to that of [5], Theorem 3.3. Corollary 2.3: The set V of solutions of the homogeneous equation Ly = 0 on R is a complex vector space of dimension n, with a basis given by {g, g , g , . . . , g(n − 1) }. If L has real coefficients (i.e., a1 , . . . , an in (2.3) are all real numbers), then the set VR of real solutions of Ly = 0 on R is a real vector space of dimension n, with a basis given by {g, g , g , . . . , g(n − 1) }. Proof: Let y be any solution of Ly = 0 on R. Let bk = y(k) (0) (0 ≤ k ≤ n − 1), and define yh by formula (2.18), with the ck given by (2.16). Then, y and yh both satisfy (2.17), whence y = yh by Theorem 2.2. We conclude by (2.18) that every solution of Ly = 0 on R is in C ∞ (R) and can be written as a linear combination of g, g , g , . . . , g(n − 1) in a unique way. (The coefficients ck in this combination are uniquely determined by the initial data bk at x = 0 through (2.16).) In particular, the n functions g, g , g , . . . , g(n − 1) are linearly independent in V and form a basis of V. If L has real coefficients, then the solutions of Ly = 0 with real initial data at x = 0 (i.e., bk = y (k) (0) ∈ R for 0 ≤ k ≤ n − 1) are real-valued. Indeed by writing y = Re y + iIm y, one gets that Im y solves the problem (2.10), whence Im y = 0. In particular, the impulsive response g is real-valued. The set VR of real solutions of Ly = 0 on R is then a real vector space of dimension n, with a basis given by {g, g , g , . . . , g(n − 1) }. 3. Explicit formulas for the impulsive response We now consider the problem of explicitly computing g for an equation of order n. By iterating (2.9), we obtain the following formula for the impulsive response of order n as an (n − 1)-times repeated integral of exponential functions 



x

gλ1 ···λn (x) =

tn−1

λn (x−tn−1 )

0

e 

eλn−1 (tn−1 −tn−2 ) · · ·

0 t3

···

e 0

λ3 (t3 −t2 )



t2

e

λ2 (t2 −t1 ) λ1 t1

e

 dt1 dt2

0

For n = 2, we have 

x

gλ1 λ2 (x) =

eλ2 (x−t) eλ1 t dt, 0

and we obtain ⎧ λx ⎨ e 1 − e λ2 x if λ1 =  λ2 , gλ1 λ2 (x) = ⎩ λλ11 x− λ2 if λ1 = λ2 . xe



 · · · dtn−2 dtn−1 . (3.1)

12

R. Camporesi For n = 3, we have  gλ1 λ2 λ3 (x) =

x

eλ3 (x−t) gλ1 λ2 (t) dt, 0

and we easily obtain the following formulas for g = gλ1 λ2 λ3 . (1) If λ1 = λ2 = λ3 (all distinct roots), then g(x) =

e λ1 x e λ2 x e λ3 x + + . (λ1 − λ2 )(λ1 − λ3 ) (λ2 − λ1 )(λ2 − λ3 ) (λ3 − λ1 )(λ3 − λ2 )

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

(2) If λ1 = λ2 = λ3 , then g(x) =

 λ3 x  1 e − eλ1 x + (λ1 − λ3 ) x eλ1 x . 2 (λ1 − λ3 )

(3) If λ1 = λ2 = λ3 , then g(x) =

1 2 λ1 x x e . 2

Note that if λ1 ∈ R and λ2, 3 = α ± iβ with α, β ∈ R, β = 0, then one computes   1 1 λ1 x αx αx g(x) = e + (α − λ1 ) e sin βx − e cos βx . (λ1 − α)2 + β 2 β

(3.2)

Example 1: Solve the initial value problem 

y  + y  = y( π2 ) =

1 sin x y  ( π2 ) =

y  ( π2 ) = 0.

Solution. The characteristic polynomial p(λ) = λ3 + λ = λ(λ2 + 1) has the roots λ1 = 0, λ2 = i, λ3 = −i. Using (3.2), we compute the impulsive response g(x) = 1 − cos x. The forcing term 1/sin x is continuous in the interval I = (0, π )  x0 = π /2. Using (2.20) with ck = 0 ∀k, we get the following solution of the proposed problem in the interval I: 

1 − cos(x − t) dt sin t π/2  x  x  x 1 cos t = dt dt − cos x dt − sin x π/2 sin t π/2 sin t π/2   t x − cos x [log sin t]xπ/2 − (x − π/2) sin x = log tan 2 π/2 x = log tan − cos x log sin x − (x − π/2) sin x. 2

y(x) =

x

International Journal of Mathematical Education in Science and Technology

13

For generic n, we get simple formulas for g = gλ1 ···λn when the roots of p(λ) are all distinct or all equal. Proposition 3.1: (1) If λi = λj for i = j (all distinct roots), then g(x) = c1 eλ1 x + c2 eλ2 x + · · · + cn eλn x , where cj = 

1 1 =  (λ − λ ) p (λ i j) i=j j

(1 ≤ j ≤ n).

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

(2) If λ1 = λ2 = ··· = λn , then g(x) =

1 x n−1 eλ1 x . (n − 1)!

(3.3)

Proof: This is readily proved by induction on n using (2.9). For (1), it is useful to remind that the function g = gλ1 ···λn is symmetric in the interchange λi ↔λj , ∀i, j.  A method for computing g in the general case is given by the following result. Theorem 3.2: Let λ1 , λ2 , . . . , λk be the distinct roots of the characteristic polynomial (2.5), of multiplicities m1 , m2 , . . . , mk , and let k

1

j =1 (λ

− λj )mj

 k   cj,mj cj,1 cj,2 = + + ··· + λ − λj (λ − λj )2 (λ − λj )mj j =1

be the partial fraction expansion of differential operator (2.3) is given by

1 p(λ)

(3.4)

over C. Then, the impulsive response of the

g(x) = G1 (x)eλ1 x + G2 (x)eλ2 x + · · · + Gk (x)eλk x ,

(3.5)

where G1 , . . . , Gk are the polynomials given by Gj (x) = cj,1 + cj,2 x + cj,3

x2 x mj −1 + · · · + cj,mj 2! (mj − 1)!

(1 ≤ j ≤ k).

Proof: This theorem can be proved by induction, though it is most easily proved using the Laplace transform or distribution theory. Indeed, one can show that the Laplace transform L of the impulsive response g is precisely the reciprocal of the characteristic polynomial: (Lg)(λ) =

1 . p(λ)

Decomposing 1/p(λ) according to (3.4) and taking the inverse Laplace transform gives the result (see [1], Equation (21) p. 81). See also [2], p. 141–142, for another proof using  distribution theory in the convolution algebra D+ . 

14

R. Camporesi

We can now use the impulsive response to compute the general solution of the homogeneous equation Ly = 0 in terms of the complex exponentials eλj x . Theorem 3.3: Every solution of Ly = 0 can be written in the form yh (x) = P1 (x)eλ1 x + P2 (x)eλ2 x + · · · + Pk (x)eλk x ,

(3.6)

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

where Pj (1 ≤ j ≤ k) is a complex polynomial of degree ≤ mj − 1. Conversely, any function of this form is a solution of Ly = 0. Proof: Let yh be a solution of Ly = 0. By Corollary 2.3, yh can be written as a linear combination of g, g , . . . , g(n − 1) , i.e., it is given by (2.18) for some complex numbers c0 , c1 , . . . , cn − 1 . Substituting (3.5) in (2.18), we get an expression of the form (3.6), where the degrees of the polynomials Pj can be less than or equal to mj − 1 (depending on the coefficients c0 , c1 , . . . , cn − 1 ). For the converse, observe that since λ1 , λ2 , . . . , λk are the distinct roots of p(λ), of multiplicities m1 , m2 , . . . , mk , p(λ) can be factored according to the formula p(λ) = (λ − λ1 )m1 (λ − λ2 )m2 · · · (λ − λk )mk . Similarly, the differential operator L factors as  L=

d − λ1 dx

m1 

d − λ2 dx

m2

 ···

d − λk dx

mk ,

(3.7)

where the order of composition of the various factors is irrelevant. We now observe that if λ ∈ C and m ∈ N+ , then 

d −λ dx

m y=0



y(x) = P (x)eλx ,

(3.8)

where P is a complex polynomial of degree ≤ m − 1. Indeed, given a function y and letting h(x) = y(x)e−λx , we compute   d − λ y (x) e−λx , h (x) = y (x)e − λy(x)e = dx  2    −λx d   2 − λ y (x) e−λx , = h (x) = y (x) − 2λy (x) + λ y(x) e dx 



 h(m) (x) =

−λx



−λx

.. . m  d − λ y (x) e−λx . dx

It follows that  m d −λ y =0 dx



h(m) = 0



h is a polynomial of degree ≤ m − 1,

as claimed. Note that the implication from left to right also follows by substituting (3.3) in (2.18), as already seen in the first part of the proof for the case of k distinct roots. By writing

International Journal of Mathematical Education in Science and Technology

15

now L in the form (3.7) and using the commutativity of the various factors and (3.8), we  see that any function of the form (3.6) satisfies Lyh = 0. Theorem 3.3 identifies another basis of the vector space V of solutions of Ly = 0 on R, namely the well-known basis given explicitly by the n functions e λj x , x e λj x , x 2 e λj x , . . . ,

x mj −1 eλj x (1 ≤ j ≤ k).

Indeed, these functions belong to V (by the converse part of the theorem), and every element of V is a linear combination of them (by (3.6)). Since V has dimension n (by Corollary 2.3), these functions must be linearly independent. When L has real coefficients, we get the well-known basis of VR given by the n functions

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

eλj x , x eλj x , . . . , x nj −1 eλj x (1 ≤ j ≤ p), eαj x cos βj x, xeαj x cos βj x, . . . , x mj −1 eαj x cos βj x (1 ≤ j ≤ q), eαj x sin βj x, xeαj x sin βj x, . . . , x mj −1 eαj x sin βj x (1 ≤ j ≤ q), where λj (1 ≤ j ≤ p) is a real root of the characteristic polynomial p(λ) of multiplicity nj , and α j ± iβ j (1 ≤ j ≤ q) is a pair of complex conjugate roots of p(λ) both of multiplicity mj , with n1 + ··· + np + 2m1 + ··· + 2mq = n. Example 2: Find the general real solution of the differential equation y (4) − 2y  + y =

1 ch3 x

.

(3.9)

Solution. The characteristic polynomial p(λ) = λ4 − 2λ2 + 1 = (λ2 − 1)2 has the roots λ1 = 1, λ2 = −1, both of multiplicity m = 2. The general solution of the homogeneous equation is then yh (x) = (ax + b) ex + (cx + d) e−x

(a, b, c, d ∈ R).

To this, we must add any particular solution of (3.9); for example, the solution yp of the initial value problem 

y (4) − 2y  + y =

1

ch3 x y(0) = y  (0) = y  (0) = y  (0) = 0. From (A-3), (A-4) and (A-5), we find the impulsive response g(x) = (− 14 + 14 x)ex + ( 14 + 14 x)e−x =

1 2

x ch x − 12 sh x.

Alternatively, use the partial fraction expansion of 1/p(λ) a c 1 b d = + , + + (λ2 − 1)2 λ − 1 (λ − 1)2 λ + 1 (λ + 1)2 compute a = −1/4, b = c = d = 1/4, and use Theorem 3.2 to get the same result.

16

R. Camporesi The forcing term

is continuous on R. By (2.8), we get ∀x ∈ R

1

ch3 x

yp (x) =

1 2



x

(x − t) ch(x − t)

0

3

ch t

dt −

1 2



x

sh(x − t) ch3 t

0

dt.

Now, use the formulas ch(x − t) = ch x ch t − sh x sh t, sh(x − t) = sh x ch t − ch x sh t. The following integrals are immediate: 

1 2

ch t

 dt = th t + c,

sh t 3

ch t

dt = −

1 2ch2 t

+ c.

Integrating by parts, we compute

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

 t

1 ch2 t

 dt = t th t − log(ch t) + c,

t

sh t ch3 t

dt = −

t 2ch2 t

+

1 th t + c. 2

Thus, yp (x) =

1 1 ch x log(ch x) − x sh x. 2 4

The final result can be written as y(x) = yh (x) +

1 ch x log(ch x), 2

by noting that the term − 14 x sh x is a solution of the homogeneous equation.

Disclosure statement No potential conflict of interest was reported by the author.

References [1] Doetsch G. Introduction to the theory and application of the Laplace transformation. New York (NY): Springer-Verlag; 1974. [2] Schwartz L. Methodes mathematiques pour les sciences physiques [Mathematical methods for physical sciences]. Enseignement des Sciences. Paris: Hermann; 1961. [3] Coddington EA, Levinson N. Theory of ordinary differential equations. New York (NY): McGraw-Hill Book Company, Inc.; 1955. [4] Coddington EA. An introduction to ordinary differential equations. Prentice-Hall Mathematics Series. Englewood Cliffs (NJ): Prentice-Hall, Inc.; 1961. [5] Camporesi R. Linear ordinary differential equations with constant coefficients. Revisiting the impulsive response method using factorization. Internat J Math Ed Sci Tech. 2011;42:497–514. [6] Robin WA. Solution of differential equations by quadratures: Cauchy integral method. Internat J Math Ed Sci Tech. 1994;25:779–790.

Appendix. An explicit formula for the impulsive response In this appendix, we give an explicit formula for the impulsive response g in the general case of k distinct roots λ1 , λ2 , . . . , λk of respective multiplicities m1 , m2 , . . . , mk , with m1 + m2 + ··· + mk = n. The following result is independent of Theorem 3.2 and can be proved without using that theorem.

International Journal of Mathematical Education in Science and Technology

17

Theorem A.1: There exist polynomials G1 , G2 , . . . , Gk , of degrees m1 − 1, m2 − 1, . . . , mk − 1, respectively, such that g(x) = G1 (x)eλ1 x + G2 (x)eλ2 x + · · · + Gk (x)eλk x .

(A-1)

More precisely, G1 is the polynomial of degree m1 − 1 given, for k ≥ 2, by m1 −1 r1

G1 (x) =

r2  

rk−2  (−1)m1 −1−rk−1 rk−1 x rk−1 ! rk−1 =0 m1 +m2 −2−r1 m3 −1+r1 −r2 m4 −1+r2 −r3

···

r1 =0 r2 =0 r3 =0

×

(λ1 −

m2 −1 λ2 )m1 +m2 −r1 −1 (λ1







k−2 −rk−1 · · · mk −1+r m3 −1 m4 −1 mk −1 λ3 )m3 +r1 −r2 (λ1 − λ4 )m4 +r2 −r3 · · · (λ1 − λk )mk +rk−2 −rk−1

.

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

(A-2) To get Gj , j ≥ 2, just exchange λ1 ↔λj and m1 ↔mj in (A-2). The proof of this result will be given elsewhere. However, the main ideas of the proof using our previous theory are indicated below for the case k = 2. Let us make explicit the formula for g given by Theorem A.1 in the cases k = 2, 3. (For k = 1, x m1 −1 , follows from formula (A-2) is somewhat degenerate. The result g = G1 eλ1 , with G1 (x) = (m 1 −1)! (3.3).) For k = 2, i.e., for two distinct roots λ1 , λ2 of multiplicities m1 , m2 , we get from (A-1)–(A-2) the impulsive response g(x) = G1 (x)eλ1 x + G2 (x)eλ2 x ,

(A-3)

where  (−1)m1 −1−r m1 + m2 − r − 2 xr , r! (λ1 − λ2 )m1 +m2 −r−1 m2 − 1 r=0

(A-4)

 (−1)m2 −1−r m1 + m2 − r − 2 xr . r! (λ2 − λ1 )m1 +m2 −r−1 m1 − 1 r=0

(A-5)

m1 −1

G1 (x) =

m2 −1

G2 (x) =

Note that G1 , G2 have degrees m1 − 1, m2 − 1, respectively, and that G2 is obtained from G1 by interchanging λ1 ↔λ2 and m1 ↔m2 , so that g is symmetric under this interchange. The simplest way to prove this is to go back to formula (3.1). Suppose we label our roots as {λ1 , . . . , λn } = {λ1 , . . . , λ1 , λ2 , . . . , λ2 }, where λ1 (resp. λ2 ) is repeated m1 (resp. m2 ) times, with m1 + m2 = n. By interchanging the order of integration in some of the nested integrals in (3.1), we arrive easily at the following formula: 

x

gλ2 ···λ2 (x − t)gλ1 ···λ1 (t) dt

g(x) = 0 x = 0

1 1 (x − t)m2 −1 eλ2 (x−t) t m1 −1 eλ1 t dt. (m2 − 1)! (m1 − 1)!

(A-6)

At this point, we use the following result, that can be proved by elementary calculus. Lemma A.2: Let p, q ∈ N and b, c ∈ C, b = c. Then, ∀x ∈ R, we have 1 bx e p!q!



x

(x − t)p t q e(c−b)t dt = P (x, p, q, b−c) ebx + P (x, q, p, c−b) ecx , 0

(A-7)

18

R. Camporesi

where P(x, p, q, a), a ∈ C \ {0}, is the polynomial of degree p in x given by P (x, p, q, a) =

  p  (−1)p−r p + q − r r=0

r!

q

xr a p+q+1−r

.

(A-8)

Using (A-7)–(A-8) in (A-6) with p = m2 − 1, q = m1 − 1 gives immediately the result (A-3)– (A-5). For k ≥ 3, we can reason in a similar way to prove Theorem A.1 by formula (3.1) and Lemma A.2. For example, for k = 3 we get g = G1 eλ1 + G2 eλ2 + G3 eλ3 , where m1 +m2 −2−r m3 −1+r−s r  (−1)m1 −1−s m2 −1 m3 −1 xs G1 (x) = s! (λ1 − λ2 )m1 +m2 −r−1 (λ1 − λ3 )m3 +r−s r=0 s=0 m1 +m2 −2−r m3 −1+r−s m1 −1 m1 −1   (−1)m1 −1−s m2 −1 m3 −1 = xs , m1 +m2 −r−1 (λ − λ )m3 +r−s s! (λ − λ ) 1 2 1 3 r=s s=0

Downloaded by [Roberto Camporesi] at 10:08 19 June 2015

m1 −1

and analogous formulae for G2 and G3 .