18.Appendix B.Ordinary differential equations

This version: 22/10/2004

Appendix B

Ordinary di erential equations

In this appendix we provide a very quick review of di erential equations at the level prerequisite to use this book. The objective, as in the other appendices, is not to provide a complete overview, but to summarise the main points. The text itself covers various aspects of di erential equations theory beyond what we say here; indeed, it is possible to think of control theory as a sub-discipline of di erential equations, although I do not like to do so. We deal with scalar equations in Section B.1, and with systems of equations, using the matrix exponential, in Section B.2.

To do this, one seeks solutions of the form eλt for λ C. Substitution into the homogeneous equations then gives the polynomial

λn + pn−1λn−1 + · · · + p1λ + p0 = 0.

This is the characteristic equation and it will have n roots, counting multiplicities. Thus it will have roots λ1, . . . , λ`, with respective multiplicities k1, . . . , k`, with k1 + · · · + k` = n. Corresponding to a root λa, one constructs ka linearly independent solutions if the root is real, and 2ka linearly independent solutions if the root is complex (a complex root and its conjugate yield the same 2ka solutions). These solutions are constructed as follows.

B.1 Construction of solutions corresponding to a root Let λ C be a root of the characteristic equation of multiplicity k.

1. If λ R then k linearly independent solutions are

Applying the construction to all roots of the characteristic equation yields n linearly independent solutions y1(t), . . . , yn(t) to the homogeneous equation. Furthermore, every solution of the homogeneous equation is a linear combination of these n linearly independent solutions. Thus we take yh to be a general linear combination of the n linearly independent solutions, with at the moment unspecified coe cients:

yh(t) = c1y1(t) + · · · + cnyn(t).

Next one determines what is called a particular solution, which we shall denote yp. A particular solution is, by definition, any solution to the di erential equation (B.1). Of course, one cannot expect to arrive at such a solution in a useful way for arbitrary righthand side functions u. What we shall describe here is an essentially ad hoc procedure, useful for certain types of functions u. This procedure typically goes under the name of the method of undetermined coe cients. The idea is that one uses the fact that certain types of functions of time—polynomial functions, trigonometric functions, and exponential functions—have their form unchanged by the act of di erentiation. Thus if one takes yp to be of this form, then

will also have the same form. Therefore, one hopes to be able to determine yp by “comparing” u to the expression (B.2). This hope is founded, and leads to the following procedure.

B.2 Method of undetermined coe cients Suppose that u(t) = t`eσt(a1 sin ωt + a2 cos ωt). Let k be the smallest integer with the property that tkytest(t) is a not solution to the homogeneous equation, where ytest is any function of t of the form

ytest(t) = P (t)eσt(α1 sin ωt + α2 cos ωt),

where α1 and α2 are arbitrary constants, and P is an arbitrary nonzero polynomial of degree at most `. Often we will have k = 0. Then seek a particular solution of the form

yp(t) = tkQ(t)eσt(b1 sin ωt + b2 cos ωt),

where Q is a polynomial of degree `. Substitute into (B.1) to determine b1, b2, and the coe cients of Q. If u is a linear combination,

u(t) = µ1u1(t) + · · · + µmum(t),

of terms of the form t`eσt(a1 sin ωt + a2 cos ωt), then the above procedure may be applied separately to each term in the linear combination separately, yielding separate particular solutions yp,1, . . . , yp,m. Then the particular solution is

Once one has a particular solution, the general solution is then simply the sum of the homogeneous and particular solution:

y(t) = yh(t) + yp(t) = c1y1(t) + · · · + cnyn(t) + yp(t).

The constants c1, . . . , cn can be determined if we have an initial value problem where the initial values y(0), y(1)(0), . . . , y(n−1)(0) are specified.

It is worth demonstrating the procedure to this point on a simple example.

B.3 Example We consider the di erential equation

... −t y(t) + 3¨y(t) + 3y˙ + y(t) = te .

First let us find a solution to the homogeneous equation

...

yh(t) + 3¨yh(t) + 3y˙h + yh(t) = 0.

The characteristic equation is λ3 + 3λ2 + 3λ + 1 = 0, which has the single root λ = −1 of multiplicity 3. Our recipe for the homogeneous solution then gives

yh(t) = c1e−t + c2te−t + c3t2e−t.

Now let us obtain a particular solution. By the above procedure, we take

ytest = (β1t + β0)e−t.

We calculate that ytest(t) and tytest identically solve the di erential equation for any β1 and β0. Also, we compute that t2ytest(t) satisfies the homogeneous provided that β1 = 0. In exasperation, we compute that the only way that t3ytest(t) will satisfy the homogeneous equation if and only if β1 = β0 = 0. Thus we have k = 3, using the notation of our procedure for obtaining particular solutions. Therefore, we seek a particular solution of the form

yp(t) = t3(q1t + q0)e−t.

Substitution into the ode gives

(6q0 + 24q1t)e−t = te−t,

from which we ascertain that q0 = 0 and q1 = 241 . Thus the particular solution is

t4e−t yp(t) = 24 .

The general solution is then

y(t) = yh(t) + yp(t) = c1e−t + c2te−t + c3t2e−t + t4e−t .

24

To determine the unknown constants c1, c2, and c3, we use the initial conditions. One may determine that

B.2 Systems of ordinary di erential equations

Next we look at first-order systems of ordinary di erential equations:

where x Rn and A Rn×n. The possibility of an additional term on the right-hand side of the form f(t) is interesting, but this is handled in the text. The solution to (B.3) is given to us by the matrix exponential. For any n × n matrix A define

eA = In + A + A2 + · · · + Ak + · · ·

2! k!

One may verify that this series converges for every A, and it is a straightforward calculation to verify that x(t) = eAtx0 is the solution to the di erential equation (B.3) with the initial condition x(0) = x0. Let us present a way to compute the matrix exponential. Corresponding to the linear system (B.3) will be n linearly independent vector solutions x1(t), . . . , xn(t). If A is diagonalisable, it is relatively easy to compute these solutions, but if A is not diagonalisable, there is a bit more work involved. Nonetheless, there is a recipe for doing this, and we present it here. Note that this recipe will always work, and it can be applied, even when A is diagonalisable—in this case it just simplifies.

First one computes the eigenvalues for A. There will be n of these in total, counting algebraic multiplicities and complex conjugate pairs. One treats each eigenvalue separately. For a real eigenvalue λ0 with algebraic multiplicity k = ma(λ0), one must compute k linearly independent solutions. For a complex eigenvalue λ0 with algebraic multiplicity ` = ma(λ0),

¯

one must compute 2` linearly independent solutions, since λ0 is also necessarily an eigenvalue with algebraic multiplicity `.

We first look at how to deal with real eigenvalues. Let λ0 be one such object with algebraic multiplicity k. It is a fact that the matrix (A − λ0In)k will have rank n − k, and so will have a kernel of dimension k by the Rank-Nullity Theorem. Let u1, . . . , uk be a basis for ker((A − λ0In)k). We call each of these vectors a generalised eigenvector. If the geometric multiplicity of λ0 is also k, then the generalised eigenvectors will simply be the usual eigenvectors. If mg(λ0) < ma(λ0) then a generalised eigenvector may or may not be an eigenvector. Corresponding to each generalised eigenvector ui, i = 1, . . . , k, we will define a

Note that because ui is a generalised eigenvector, the infinite series exp((A − λ0In)t)ui will have only a finite number of terms—at most k in fact. Indeed we have

since the remaining terms in the series will be zero. In any case, it turns out that the k vector functions x1(t), . . . , xk(t) so constructed will be linearly independent solutions of (B.3). This tells us how to manage the real case.

The characteristic polynomial for A is PA(λ) = λ2 + 4λ + 3 = (λ + 1)(λ + 3). Thus the eigenvalues for A are λ1 = −1 and λ2 = −3. Since the eigenvalues are distinct, the algebraic and geometric multiplicities will be equal, and the generalised eigenvectors will simply be eigenvectors. An eigenvector for λ1 = −1 is u1 = (2, 3) and an eigenvector for λ2 = −3 is u2 = (1, 1). Our recipe then gives two linearly independent solutions as

2. A more interesting case is the following:

00 −1

Since the matrix is upper triangular, the eigenvalues are the diagonal elements: λ1 = −2 and λ2 = −1. The algebraic multiplicity of λ1 is 2. However, we readily see that dim(ker(A − λ1I3)) = 1 and so the geometric multiplicity is 1. So we need to compute generalised eigenvectors in this case. We have

and the generalised eigenvectors span the kernel of this matrix, and so we may take u1 = (1, 0, 0) and u2 = (0, 1, 0) as generalised eigenvectors. Applying the formula (B.4) gives

e−2t

=0 0

and

0

Finally we determine that u3 = (0, 0, 1) is an eigenvector corresponding to λ2 = −1, and so this gives the solution

0

x3(t) = 0 . e−t

Now let us look at the complex case. Thus let λ0 be a complex eigenvalue (with nonzero

¯

imaginary part) of algebraic multiplicity `. This means that λ0 will also be an eigenvalue of algebraic multiplicity ` since A, and hence PA(λ), is real. Thus we need to find 2` linearly independent solutions. We do this by following the exact same idea as in the real case, except that we think of A as being a complex matrix for the moment. In this case it is still true that the matrix (A − λ0In)` will have an `-dimensional kernel, and we can take vectors u1, . . . , u` as a basis for this kernel. Note, however, that since (A−λ0In)` is complex, these vectors will also be complex. But the procedure is otherwise identical to the real case. One then constructs ` complex vector functions

Each such complex vector function will be a sum of its real and imaginary parts: zj (t) =

xj (t) + iyj (t). It turns out that the 2` real vector functions x1(t), . . . , x`(t), y1(t), . . . , y`(t) are linearly independent solutions to (B.3).

Let’s see how this works in some examples.

B.5 Examples 1. An example with complex roots is

00 −2

The characteristic polynomial is PA(λ) = r3 + 4r2 + 6r + 4. One ascertains that the

from which we glean that an eigenvector is u1 = (−i, 1, 0). Using (B.5) the complex solution is then

Corresponding to the real eigenvalue λ3 we readily determine that

0

x2 = e−2t 0

1

is a corresponding solution. This gives three linearly independent real solutions x1(t), y1(t), and x2(t).

2.Let us look at a complex root with nontrivial multiplicity. The smallest matrix to possess such a feature will be one which is 4 × 4, and about the simplest example for which the geometric multiplicity is less than the algebraic multiplicity is

00 −1 0

The eigenvalues are λ1 = i and λ2 = −i, both with algebraic multiplicity 2. One readily determines that the kernel of A−iI4 is one-dimensional, and so the geometric multiplicity of these eigenvalues is just 1. Thus we need to compute complex generalised eigenvectors. I have used Mathematica® for the computations below. We compute

and one checks that u1 = (0, 0, −i, 1) and u2 = (−i, 1, 0, 0) are two linearly independent generalised eigenvectors. We compute

We now determine the two linearly independent real solutions corresponding to u1. We have

z1(t) = eit(u1 + t(A − iI4)u1)

0−i

=eit 0 + teit 1

−i 0

10

the resulting matrix is an example of a fundamental matrix. Generally, a fundamental matrix is any n×n matrix function of t whose columns form n linearly independent solutions to (B.3). What we have done above is give a recipe for computing a fundamental matrix (there are an infinite number of these). The following result connects the construction of a fundamental matrix with the matrix exponential.

B.6 Theorem Given any fundamental matrix X(t) we have eAt = X(t)X−1(0).

Thus, once we have a fundamental matrix, the computation of the matrix exponential is just algebra, although computing inverses of matrices of any size is a task best left to the computer.

Let’s work this out for our four examples.

B.7 Examples 1. If

1Note that the solutions x1, . . . , xn are those obtained from both real and complex eigenvalues. Therefore, the solutions denoted above as “yi(t)” for complex eigenvalues will be included in the n linearly independent solutions, except now I am calling everything xj (t).