3 Systems of ODEs

[Go to the Notes on Diffy Qs home page] [PDF version]

Chapter 3
Systems of ODEs

3.1 Introduction to systems of ODEs

Note: 1 lecture, §4.1 in EP

Often we do not have just one dependent variable and one equation. And as we will see, we may end up with systems of several equations and several dependent variables even if we start with a single equation.

If we have several dependent variables, suppose $y 1$ , $y 2$ , …, $y n$ we can have a differential equation involving all of them and their derivatives. For example, $′′ ′ ′ y1 = f(y1,y2,y1,y2,x)$ . Usually, when we have two dependent variables we would have two equations such as

′′ ′ ′ y1 = f1(y1,y2,y1,y2,x), y′′2 = f2(y′1,y′2,y1,y2,x),

for some functions $f1$ and $f2$ . We call the above a system of differential equations. More precisely, it is a second order system. Sometimes a system is easy to solve by solving for one variable and then for the second variable.

Example 3.1.1: Take the first order system

′ y1 = y1, y′2 = y1 - y2,

with initial conditions of the form $y1(0) = 1$ , $y2(0 ) = 2$ .

We note that $y = C ex 1 1$ is the general solution of the first equation. We can then plug this $y 1$ into the second equation and get the equation $′ x y2 = C 1e - y2$ , which is a linear first order equation that is easily solved for $y2$ . By the method of integrating factor we get

exy = C-1e2x + C , 2 2 2

or $y2 = C21 ex + C2e-x$ . The general solution to the system is, therefore,

y1 = C1ex, C1 y2 = --ex + C2e-x. 2

We can now solve for $C1$ and $C2$ given the initial conditions. We substitute $x = 0$ and find that $C1 = 1$ and $C 2 = 3∕2$ .

Generally, we will not be so lucky to be able to solve like in the first example, and we will have to solve for all variables at once.

As an example application, let us think of mass and spring systems again. Suppose we have one spring with constant $k$ but two masses $m 1$ and $m 2$ . We can think of the masses as carts, and we will suppose that they ride along with no friction. Let $x1$ be the displacement of the first cart and $x2$ be the displacement of the second cart. That is, we put the two carts somewhere with no tension on the spring, and we mark the position of the first and second cart and call those the zero position. That is, $x1 = 0$ is a different position on the floor than the position corresponding to $x2 = 0$ . The force exerted by the spring on the first cart is $k(x 2 - x1)$ , since $x2 - x1$ is how far the string is stretched (or compressed) from the rest position. The force exerted on the second cart is the opposite, thus the same thing with a negative sign. Newton’s second law states that force equals mass times acceleration.

m 1x′′ = k(x2 - x1), 1′′ m 2x2 = -k(x2 - x1).

In this system we cannot solve for the $x 1$ variable separately. That we must solve for both $x 1$ and $x 2$ at once is intuitively obvious, since where the first cart goes depends exactly on where the second cart goes and vice-versa.

Before we talk about how to handle systems, let us note that in some sense we need only consider first order systems. Take an $th n$ order differential equation

y(n) = F (y(n-1),...,y′,y,x).

Define new variables $u ,...,u 1 n$ and write the system

u′1 = u2 ′ u2 = u3 .. . u′n-1 = un ′ un = F(un,un-1,...,u 2,u 1,x ).

Now try to solve this system for $u1$ , $u 2$ , …, $un$ . Once you have solved for the $u$ ’s, you can discard $u 2$ through $un$ and let $y = u1$ . We note that this $y$ solves the original equation.

A similar process can be followed for a system of higher order differential equations. For example, a system of $k$ differential equations in $k$ unknowns, all of order $n$ , can be transformed into a first order system of $n × k$ equations and $n × k$ unknowns.

Example 3.1.2: Sometimes we can use this idea in reverse as well. Let us take the system

′ ′ x = 2y - x, y = x,

where the independent variable is $t$ . We wish to solve for the initial conditions $x(0) = 1$ , $y(0) = 0$ .

We first notice that if we differentiate the first equation once we get $′′ ′ y = x$ and now we know what $′ x$ is in terms of $x$ and $y$ .

y′′ = x′ = 2y - x = 2y - y′.

So we now have an equation $y′′ + y ′ - 2y = 0$ . We know how to solve this equation and we find that $y = C e-2t + C et 1 2$ . Once we have $y$ we can plug in to get $x$ .

x = y′ = -2C 1e-2t + C2et.

We solve for the initial conditions $1 = x(0 ) = - 2C1 + C2$ and $0 = y(0) = C1 + C2$ . Hence, $C 1 = - C2$ and $1 = 3C 2$ . So $-1 C1 = ∕3$ and $1 C2 = ∕3$ . Our solution is

2e-2t + et -e-2t + et x = --------, y = ---------. 3 3

Exercise 3.1.1: Plug in and check that this really is the solution.

It is useful to go back and forth between systems and higher order equations for other reasons. For example, the ODE approximation methods are generally only given as solutions for first order systems. It is not very hard to adapt the code for the Euler method for a first order equation to first order systems. We essentially just treat the dependent variable not as a number but as a vector. In many mathematical computer languages there is almost no distinction in syntax.

In fact, this is what IODE was doing when you had it solve a second order equation numerically in the IODE Project III if you have done that project.

The above example was what we will call a linear first order system, as none of the dependent variables appear in any functions or with any higher powers than one. It is also autonomous as the equations do not depend on the independent variable $t$ .

For autonomous systems we can easily draw the so-called direction field or vector field. That is, a plot similar to a slope field, but instead of giving a slope at each point, we give a direction (and a magnitude). The previous example $′ x = 2y - x$ , $′ y = x$ says that at the point $(x,y)$ the direction in which we should travel to satisfy the equations should be the direction of the vector $(2y - x,x)$ with the speed equal to the magnitude of this vector. So we draw the vector $(2y - x,x)$ based at the point $(x,y)$ and we do this for many points on the $xy$ -plane. We may want to scale down the size of our vectors to fit many of them on the same direction field. See Figure 3.1.

We can now draw a path of the solution in the plane. That is, suppose the solution is given by $x = f(t)$ , $y = g(t)$ , then we can pick an interval of $t$ (say $0 ≤ t ≤ 2$ for our example) and plot all the points $(f(t),g(t))$ for $t$ in the selected range. The resulting picture is usually called the phase portrait (or phase plane portrait). The particular curve obtained we call the trajectory or solution curve. An example plot is given in Figure 3.2. In this figure the line starts at $(1,0 )$ and travels along the vector field for a distance of 2 units of $t$ . Since we solved this system precisely we can compute $x(2)$ and $y(2)$ . We get that $x(2) ≈ 2.4 75$ and $y(2) ≈ 2.457$ . This point corresponds to the top right end of the plotted solution curve in the figure.

Figure 3.1: The direction field for $x′ = 2y - x$ , $y′ = x$ .

Figure 3.2: The direction field for $x′ = 2y - x$ , $y′ = x$ with the trajectory of the solution starting at $(1,0)$ for $0 ≤ t ≤ 2$ .

Notice the similarity to the diagrams we drew for autonomous systems in one dimension. But now note how much more complicated things become if we allow just one more dimension.

Also note that we can draw phase portraits and trajectories in the $xy$ -plane even if the system is not autonomous. In this case however we cannot draw the direction field, since the field changes as $t$ changes. For each $t$ we would get a different direction field.

3.1.1 Exercises

Exercise 3.1.2: Find the general solution of $x′1 = x2 - x1 + t$ , $x′2 = x2$ .

Exercise 3.1.3: Find the general solution of $′ t x1 = 3x1 - x2 + e$ , $′ x2 = x1$ .

Exercise 3.1.4: Write $′′ ′ ay + by + cy = f(x)$ as a first order system of ODEs.

Exercise 3.1.5: Write $x′′ + y2y′ - x3 = sin(t)$ , $y′′ + (x ′ + y′)2 - x = 0$ as a first order system of ODEs.

3.2 Matrices and linear systems

Note: 1 and a half lectures, first part of §5.1 in EP

3.2.1 Matrices and vectors

Before we can start talking about linear systems of ODEs, we will need to talk about matrices, so let us review these briefly. A matrix is an $m × n$ array of numbers ( $m$ rows and $n$ columns). For example, we denote a $3 × 5$ matrix as follows

⌊ ⌋ |||||a11 a12 a13 a14 a15||||| A = |||⌈a21 a22 a23 a24 a25|||⌉. a31 a32 a33 a34 a35

By a vector we will usually mean a column vector which is an $n × 1$ matrix. If we mean a row vector we will explicitly say so (a row vector is a $1 × n$ matrix). We will usually denote matrices by upper case letters and vectors by lower case letters with an arrow such as $⃗x$ or $⃗b$ . By $⃗0$ we will mean the vector of all zeros.

It is easy to define some operations on matrices. Note that we will want $1 × 1$ matrices to really act like numbers, so our operations will have to be compatible with this viewpoint.

First, we can multiply by a scalar (a number). This means just multiplying each entry by the same number. For example,

[ ] [ ] 2 1 2 3 = 2 4 6 . 4 5 6 8 10 12

Matrix addition is also easy. We add matrices element by element. For example,

[ ] [ ] [ ] 1 2 3 1 1 - 1 2 3 2 4 5 6 + 0 2 4 = 4 7 10 .

If the sizes do not match, then addition is not defined.

If we denote by 0 the matrix of with all zero entries, by $c$ , $d$ some scalars, and by $A$ , $B$ , $C$ some matrices, we have the following familiar rules.

A + 0 = A = 0 + A, A + B = B + A, (A + B) + C = A + (B + C ), c(A + B) = cA + cB, (c + d)A = cA + dA.

Another operation which is useful for matrices is the so-called transpose. This operation just swaps rows and columns of a matrix. The transpose of $A$ is denoted by $AT$ . Example:

⌊ ⌋ [ ]T |||1 4||| 1 2 3 = ||||2 5|||| 4 5 6 |⌈3 6|⌉

3.2.2 Matrix multiplication

Let us now define matrix multiplication. First we define the so-called dot product (or inner product) of two vectors. Usually this will be a row vector multiplied with a column vector of the same size. For the dot product we multiply each pair of entries from the first and the second vector and we sum these products. The result is a single number. For example,

⌊ ⌋ [ ] ||||b1|||| a1 a2 a3 ⋅||||b2|||| = a1b1 + a2b2 + a3b3. ⌈b3⌉

And similarly for larger (or smaller) vectors.

Armed with the dot product we can define the product of matrices. First let us denote by $rowi(A)$ the $ith$ row of $A$ and by $colu mnj(A)$ the $jth$ column of $A$ . Now for an $m × n$ matrix $A$ and an $n × p$ matrix $B$ we can define the product $AB$ . We let $AB$ be an $m × p$ matrix whose $ijth$ entry is

ro wi(A ) ⋅ colum nj(B).

Do note how the sizes match up. Example:

⌊ ⌋
[1 2 3] |||1 0 -1|||
|||||1 1 1 ||||| =
4 5 6 ⌈1 0 0 ⌉
[ ] [ ]
= 1 ⋅ 1 + 2 ⋅ 1 + 3 ⋅ 1 1 ⋅ 0 + 2 ⋅ 1 + 3 ⋅ 0 1 ⋅ (-1 ) + 2 ⋅ 1 + 3 ⋅ 0 = 6 2 1
4 ⋅ 1 + 5 ⋅ 1 + 6 ⋅ 1 4 ⋅ 0 + 5 ⋅ 1 + 6 ⋅ 0 4 ⋅ (-1 ) + 5 ⋅ 1 + 6 ⋅ 0 15 5 1

For multiplication we will want an analogue of a 1. This is the so-called identity matrix. The identity matrix is a square matrix with 1s on the main diagonal and zeros everywhere else. It is usually denoted by $I$ . For each size we have a different identity matrix and so sometimes we may denote the size as a subscript. For example, the $I3$ would be the $3 × 3$ identity matrix

⌊|1 0 0⌋| |||| |||| I = I3 = |||⌈0 1 0|||⌉. 0 0 1

We have the following rules for matrix multiplication. Suppose that $A$ , $B$ , $C$ are matrices of the correct sizes so that the following make sense. Let $α$ denote a scalar (number).

A(BC ) = (AB )C, A(B + C ) = AB + AC, (B + C)A = BA + CA, α(AB ) = (αA )B = A(αB ), IA = A = AI.

A few warnings are in order however.

(i): $AB ⇔ BA$ in general (it may be true by fluke sometimes). That is, matrices do not commute.
(ii): $AB = AC$ does not necessarily imply $B = C$ even if $A$ is not 0.
(iii): $AB = 0$ does not necessarily mean that $A = 0$ or $B = 0$ .

For the last two items to hold we would need to essentially “divide” by a matrix. This is where matrix inverse comes in. Suppose that $A$ is an $n × n$ matrix and that there exists another $n × n$ matrix $B$ such that

AB = I = BA.

Then we call $B$ the inverse of $A$ and we denote $B$ by $A -1$ . If the inverse of $A$ exists, then we call $A$ invertible. If $A$ is not invertible we say $A$ is singular or just say it is not invertible.

If $A$ is invertible, then $AB = AC$ does imply that $B = C$ (this also means that the inverse is unique). We just multiply both sides by $A-1$ to get $A -1AB = A -1AC$ or $IB = IC$ or $B = C$ . It is also not hard to see that $(A-1)-1 = A$ .

3.2.3 The determinant

We can now talk about determinants of square matrices. We define the determinant of a $1 × 1$ matrix as the value of its own entry. For a $2 × 2$ matrix we define

([ ]) det a b = ad - bc. c d

Before trying to compute determinant for larger matrices, let us first note the meaning of the determinant. Consider an $n × n$ matrix as a mapping of $n ℝ$ to $n ℝ$ . For example, a $2 × 2$ matrix $A$ is a mapping of the plane where $⃗x$ gets sent to $A⃗x$ . Then the determinant of $A$ is then the factor by which the volume of objects gets changed. For example, if we take the unit square (square of sides 1) in the plane, then $A$ takes the square to a parallelogram of area $|det(A )|$ . The sign of $det(A )$ denotes changing of orientation (if the axes got flipped). For example,

[ ] A = 1 1 . -1 1

Then $det(A ) = 1 + 1 = 2$ . Now let us see where the square with vertices $(0,0)$ , $(1,0 )$ , $(0,1 )$ , and $(1,1)$ gets sent. Obviously $(0,0)$ gets sent to $(0,0)$ . Now

[ 1 1 ] [1] [ 1] [ 1 1] [0] [1 ] [ 1 1] [1] [2 ] = , = , = . -1 1 0 -1 - 1 1 1 1 -1 1 1 0

So it turns out that the image of the square is another square. This one has a side of length $√ -- 2$ and is therefore of area 2.

If you think back to high school geometry, you may have seen a formula for computing the area of a parallelogram with vertices $(0,0)$ , $(a,c)$ , $(b,d)$ and $(a + b,c + d )$ . And it is precisely

|| ([ ])|| ||d et a b ||. || c d ||

The vertical lines here mean absolute value. The matrix $[a b] c d$ carries the unit square to the given parallelogram.

Now we can define the determinant for larger matrices. We define $Aij$ as the matrix $A$ with the $th i$ row and the $jth$ column deleted. To compute the determinant of a matrix, pick one row, say the $ith$ row and compute.

|------------------------------| | ∑n i+j | | det(A ) = (-1 ) aijdet(Aij). | -----------j=1-------------------

For example, for the first row we get

( ||+a det(A ) if n is odd, det(A) = a11det(A 11) - a12det(A 12) + a13det(A 13) - ⋅⋅⋅{|| 1n 1n (- a1ndet(A 1n) if n even.

We alternately add and subtract the determinants of the submatrices $Aij$ for a fixed $i$ and all $j$ . For example, for a $3 × 3$ matrix, picking the first row, we would get $det(A) = a11 det(A 11) - a12 det(A 12) + a13 det(A 13)$ . For example,

( ⌊ ⌋) || ||1 2 3 |||| ([ ]) ( [ ]) ([ ]) det|||| ||||4 5 6 |||||||| = 1 ⋅ det 5 6 - 2 ⋅ det 4 6 + 3 ⋅ d et 4 5 ||( ||⌈ ||⌉||) 8 9 7 9 7 8 7 8 9 = 1(5 ⋅ 9 - 6 ⋅ 8) - 2 (4 ⋅ 9 - 6 ⋅ 7 ) + 3(4 ⋅ 8 - 5 ⋅ 7) = 0.

The numbers $(-1 )i+jdet(Aij)$ are called cofactors of the matrix and this way of computing the determinant is called the cofactor expansion. It is also possible to compute the determinant by expanding along columns (picking a column instead of a row above).

Note that a common notation for the determinant is a pair of vertical lines.

|| || ([ ]) ||a b|| a b ||c d|| = det c d .

I personally find this notation confusing since vertical lines for me usually mean a positive quantity, while determinants can be negative. So I will not ever use this notation in these notes.

One of the most important properties of determinants (in the context of this course) is the following theorem.

Theorem 3.2.1. An $n × n$ matrix $A$ is invertible if and only if $det(A) ⇔ 0$ .

In fact, we have a formula for the inverse of a $2 × 2$ matrix

[ ]-1 [ ] a b = --1---- d - b . c d ad - bc -c a

Notice the determinant of the matrix in the denominator of the fraction. The formula only works if the determinant is nonzero, otherwise we are dividing by zero.

3.2.4 Solving linear systems

One application of matrices we will need is to solve systems of linear equations. This may be best shown by example. Suppose that we have the following system of linear equations

2x1 + 2x2 + 2x3 = 2, x1 + x2 + 3x3 = 5, x1 + 4x2 + x3 = 10.

Without changing the solution, we could swap equations in this system, we could multiply any of the equations by a nonzero number, and we could add a multiple of one equation to another equation. It turns out these operations always suffice to find a solution.

It is easier to write the system as a matrix equation. Note that the system can be written as

⌊ ⌋ ⌊ ⌋ ⌊ ⌋ |||2 2 2||| |||x1||| |||2 ||| ||||1 1 3|||| ||||x2||||= ||||5 ||||. |⌈1 4 1|⌉ |⌈x |⌉ |⌈10 |⌉ 3

To solve the system we put the coefficient matrix (the matrix on the left hand side of the equation) together with the vector on the right and side and get the so-called augmented matrix

⌊ ⌋ ||2 2 2 2 || |||||1 1 3 5 |||||. |⌈1 4 1 10|⌉

We then apply the following three elementary operations.

(i): Swap two rows.
(ii): Multiply a row by a nonzero number.
(iii): Add a multiple of one row to another row.

We will keep doing these operations until we get into a state where it is easy to read off the answer or until we get into a contradiction indicating no solution, for example if we come up with an equation such as $0 = 1$ .

Let us work through the example. First multiply the first row by $1∕2$ to obtain

⌊|1 1 1 1 ⌋| ||||| ||||| ||⌈1 1 3 5 ||⌉. 1 4 1 10

Now subtract the first row from the second and third row.

⌊ ⌋ ||1 1 1 1|| ||||0 0 2 4|||| ||⌈ ||⌉ 0 3 0 9

Multiply the last row by $1∕3$ and the second row by $1∕2$ .

⌊ ⌋ ||||1 1 1 1|||| ||||0 0 1 2|||| ⌈0 1 0 3⌉

Swap rows 2 and 3.

⌊ ⌋ |||1 1 1 1||| ||||0 1 0 3|||| |⌈0 0 1 2|⌉

Subtract the last row from the first, then subtract the second row from the first.

⌊||1 0 0 -4⌋|| ||||0 1 0 3|||| ||⌈ ||⌉ 0 0 1 2

If we think about what equations this augmented matrix represents, we see that $x1 = -4$ , $x 2 = 3$ , and $x3 = 2$ . We try this solution in the original system and, voilà, it works!

Exercise 3.2.1: Check that this solution works.

If we write this equation in matrix notation as

A⃗x = ⃗b,

where $A$ is the matrix $[2 2 2] 11 1 4 31$ and $⃗b$ is the vector $[2] 150$ . The solution can be also computed with the inverse,

-1 -1 ⃗x = A A⃗x = A ⃗b.

One last note to make about linear systems of equations is that it is possible that the solution is not unique (or that no solution exists). It is easy to tell if a solution does not exist. If during the row reduction you come up with a row where all the entries except the last one are zero (the last entry in a row corresponds to the right hand side of the equation) the system is inconsistent and has no solution. For example if for a system of 3 equations and 3 unknowns you find a row such as $[0 0 0 1]$ in the augmented matrix, you know the system is inconsistent.

You generally try to use row operations until the following conditions are satisfied. The first nonzero entry in each row is called the leading entry.

(i): There is only one leading entry in each column.
(ii): All the entries above and below a leading entry are zero.
(iii): All leading entries are 1.

Such a matrix is said to be in reduced row echelon form. The variables corresponding to columns with no leading entries are said to be free variables. Free variables mean that we can pick those variables to be anything we want and then solve for the rest of the unknowns.

Example 3.2.1: The following augmented matrix is in reduced row echelon form.

⌊|1 2 0 3⌋| |||| |||| |||⌈0 0 1 1|||⌉ 0 0 0 0

If the variables are named $x1$ , $x2$ , and $x3$ , then $x2$ is the free variable and $x1 = 3 - 2x 2$ and $x3 = 1$ .

On the other hand if during the row reduction process you come up with the matrix

⌊ ⌋ |||1 2 1 3 3||| ||||0 0 1 1||||, |⌈0 0 0 3|⌉

there is no need to go further. The last row corresponds to the equation $0x 1 + 0x2 + 0x 3 = 3$ which is preposterous. Hence, no solution exists.

3.2.5 Computing the inverse

If the coefficient matrix is square and there exists a unique solution $⃗x$ to $A⃗x = ⃗b$ for any $⃗b$ , then $A$ is invertible. In fact by multiplying both sides by $A -1$ you can see that $⃗x = A -1⃗b$ . So it is useful to compute the inverse, if you want to solve the equation for many different right hand sides $⃗ b$ .

The $2 × 2$ inverse is basically given by a formula, but it is not hard to also compute inverses of larger matrices. While we will not have too much occasion to compute inverses for larger matrices than $2 × 2$ by hand, let us touch on how to do it. Finding the inverse of $A$ is actually just solving a bunch of linear equations. If you can solve $A ⃗x = ⃗e k k$ where $⃗e k$ is the vector with all zeros except a 1 at the $kth$ position, then the inverse is the matrix with the columns $⃗xk$ for $k = 1,...,n$ (exercise: why?). Therefore, to find the inverse we can write a larger $n × 2n$ augmented matrix $[A I]$ , where $I$ is the identity. If you do row reduction and put the matrix in reduced row echelon form, then the matrix will be of the form $[I A -1]$ if and only if $A$ is invertible, so you can just read off the inverse $-1 A$ .

3.2.6 Exercises

Exercise 3.2.2: Solve $[ ] [ ] 13 2 4 ⃗x = 56$ by using matrix inverse.

Exercise 3.2.3: Compute determinant of $[9 -2 -6] -8 3 6 10 -2 -6$ .

Exercise 3.2.4: Compute determinant of $[ ] 14 20 35 1 0 6 0 7 0 8 0 10 1$ . Hint: expand along the proper row or column to make the calculations simpler.

Exercise 3.2.5: Compute inverse of $[ ] 1 2 3 10 1 1 10$ .

Exercise 3.2.6: For which $h$ is $[1 2 3] 47 5 8 6h$ not invertible? Is there only one such $h$ ? Are there several? Infinitely many.

Exercise 3.2.7: For which $h$ is $[h 1 1] 0 h 0 1 1 h$ not invertible? Find all such $h$ .

Exercise 3.2.8: Solve $[ ] [] -98 -23 -66 ⃗x = 12 10 -2 -6 3$ .

Exercise 3.2.9: Solve $[5 3 7] [2] 86 43 4 3 ⃗x = 00$ .

Exercise 3.2.10: Solve $[ ] [ ] 33 23 3 3 0 3 20 02 23 4 4 2 3 ⃗x = 41$ .

3.3 Linear systems of ODEs

Note: less than 1 lecture, second part of §5.1 in EP

First let us talk about matrix or vector valued functions. This is essentially just a matrix whose entries depend on some variable. Let us say the independent variable is $t$ . Then a vector valued function $⃗x(t)$ is really something like

⌊ ⌋ |||x1(t)||| ||||x2(t)|||| ⃗x(t) = ||||| .. ||||| . |||⌈ . |||⌉ xn(t)

Similarly a matrix valued function is something such as

⌊ ⌋ |||a11(t) a12(t) ⋅⋅⋅ a1n(t)||| ||||a21(t) a22(t) ⋅⋅⋅ a2n(t)|||| A (t) = |||| . . . . ||||. ||||⌈ .. .. .. .. ||||⌉ an1(t) an2(t) ⋅⋅⋅ ann(t)

We can talk about the derivative $A′(t)$ or $dA- dt$ and this is just the matrix valued function whose $ijth$ entry is $′ aij(t)$ .

Rules of differentiation of matrix valued functions are similar to rules for normal functions. Let $A$ and $B$ be matrix valued functions. Let $c$ a scalar and $C$ be a constant matrix. Then

′ ′ ′ (A + B ) = A + B (AB )′ = A′B + AB ′ ′ ′ (cA ) = cA (CA )′ = CA ′ ′ ′ (AC ) = A C

Do note the order in the last two expressions.

A first order linear system of ODEs is a system which can be written as

⃗x ′(t) = P (t)⃗x(t) + ⃗f(t).

Where $P$ is a matrix valued function, and $⃗x$ and $⃗ f$ are vector valued functions. We will often suppress the dependence on $t$ and only write $′ ⃗ ⃗x = P ⃗x +f$ . A solution is of course a vector valued function $⃗x$ satisfying the equation.

For example, the equations

x′1 = 2tx1 + etx2 + t2, ′ x1 t x2 = -t - x2 + e ,

can be written as

[2t et] [t2] ⃗x′ = 1 ⃗x + t . ∕t - 1 e

We will mostly concentrate on equations that are not just linear, but are in fact constant coefficient equations. That is, the matrix $P$ will be a constant and not depend on $t$ .

When $⃗f = ⃗0$ (the zero vector), then we say the system is homogeneous. For homogeneous linear systems we still have the principle of superposition, just like for single homogeneous equations.

Theorem 3.3.1 (Superposition). Let $⃗x′ = P ⃗x$ be a linear homogeneous system of ODEs. Suppose that $⃗x1,...,⃗xn$ are $n$ solutions of the equation, then

⃗x = c1⃗x1 + c2x⃗2 + ⋅⋅⋅ + cn⃗xn,

(3.1)

is also a solution. If furthermore this is a system of $n$ equations ( $P$ is $n × n$ ), and $⃗x1,...,⃗xn$ are linearly independent. Then every solution can be written as (3.1).

Linear independence for vector valued functions is essentially the same as for normal functions. $⃗x1,...,⃗xn$ are linearly independent if and only if

c1⃗x1 + c2⃗x2 + ⋅⋅⋅ + cn⃗xn = ⃗0

has only the solution $c1 = c2 = ⋅⋅⋅ = cn = 0$ .

The linear combination $c⃗x + c ⃗x + ⋅⋅⋅ + c ⃗x 1 1 2 2 n n$ could always be written as

X(t)⃗c,

where $X (t)$ is the matrix with columns $⃗x1,...,⃗xn$ , and $⃗c$ is the column vector with entries $c1,...,cn$ . $X (t)$ is called the fundamental matrix, or fundamental matrix solution.

To solve nonhomogeneous first order linear systems. We apply the same technique as we did before.

Theorem 3.3.2. Let $′ ⃗x = P⃗x + ⃗f$ be a linear system of ODEs. Suppose $⃗xp$ is one particular solution. Then every solution can be written as

⃗x = ⃗xc + ⃗xp,

where $⃗xc$ is a solution to the associated homogeneous equation ( $′ ⃗x = P⃗x$ ).

So the procedure will be exactly the same. We find a particular solution to the nonhomogeneous equation, then we find the general solution to the associated homogeneous equation and we add the two.

Alright, suppose you have found the general solution $⃗x′ = Px⃗+ f⃗$ . Now you are given an initial condition of the form $⃗ ⃗x(t0) = b$ for some constant vector $⃗ b$ . Suppose that $X(t)$ is the fundamental matrix solution of the associated homogeneous equation (i.e. columns of $X$ are solutions). The general solution is written as

x⃗(t) = X(t)⃗c + xp(t).

Then we are seeking a vector $⃗c$ such that

⃗b = ⃗x(t0) = X (t0)⃗c + xp(t0).

Or in other words we are solving the nonhomogeneous system of linear equations

X (t0)⃗c = ⃗b - xp(t0)

for $⃗c$ .

Example 3.3.1: In §3.1 we solved the following system

′ x1 = x 1, x′2 = x 1 - x2.

with initial conditions $x1(0) = 1$ , $x2(0) = 2$ .

This is a homogeneous system, so $f⃗= ⃗0$ . We write the system as

[ ] [ ] ′ 1 0 1 ⃗x = 1 - 1 ⃗x, ⃗x(0) = 2 .

We found the general solution was $x1 = c1et$ and $x 2 = c12 et + c2e-t$ . Hence in matrix notation, the fundamental matrix solution is

[ t ] X (t) = e 0 . 12et e-t

It is not hard to see that the columns of this matrix are linearly independent. To see this, just plug in $t = 0$ and note that the two constant vectors are already linearly independent here.

Hence to solve the initial problem we solve the equation

X(0)⃗c = ⃗b,

or in other words,

[1 0] [1] 1 ⃗c = . 2 1 2

After a single elementary row operation we find that $[ ] ⃗c = 31∕2$ . Hence our solution is

[ t ] [ ] [ t ] ⃗x(t) = X (t)⃗c = e1 t 0-t 13 = 1 te 3 -t . 2e e 2 2e + 2e

This agrees with our previous solution.

3.3.1 Exercises

Exercise 3.3.1: Write the system $x′1 = 2x 1 - 3tx2 + sint$ , $x′2 = etx1 + 3x2 + cost$ as in the form $⃗x′ = P(t)⃗x + ⃗f(t)$ .

Exercise 3.3.2: a) Verify that the system $′ [1 3] ⃗x = 3 1 ⃗x$ has the two solutions $[1] 4t 1 e$ and $[1 ] -2t -1 e$ . b) Write down the general solution. c) Write down the general solution in the form $x1 =?$ , $x 2 = ?$ (i.e. write down a formula for each element of the solution).

Exercise 3.3.3: Verify that $[1]et 1$ and $[1 ]et -1$ are linearly independent. Hint: Just plug in $t = 0$ .

Exercise 3.3.4: Verify that $[1] t 10 e$ and $[1 ] t -11 e$ and $[1] 2t -11e$ are linearly independent. Hint: You must be a bit more tricky than in the previous exercise.

Exercise 3.3.5: Verify that $[t] t2$ and $[t3] t4$ are linearly independent.

3.4 Eigenvalue method

Note: 2 lectures, §5.2 in EP

In this section we will learn how to solve linear homogeneous constant coefficient systems of ODEs by the eigenvalue method. Suppose you have a linear constant coefficient homogeneous system

⃗x′ = P ⃗x.

Now suppose we try to adapt the method for single constant coefficient equations by trying the function $λt e$ . However, $⃗x$ is a vector. So we try $λt ⃗ve$ , where $⃗v$ is an arbitrary constant vector. We plug into the equation to get

λ⃗veλt = P⃗veλt.

We divide by $eλt$ and notice that we are looking for a $λ$ and $⃗v$ that satisfy the equation

λ⃗v = P⃗v.

To solve this equation we need a little bit more linear algebra which we review now.

3.4.1 Eigenvalues and eigenvectors of a matrix

Let $A$ be a square constant matrix. Suppose there is a scalar $λ$ and a nonzero vector $⃗v$ such that

A⃗v = λ⃗v.

We then call $λ$ an eigenvalue of $A$ and $⃗v$ is called the corresponding eigenvector.

Example 3.4.1: The matrix $[2 1] 0 1$ has an eigenvalue of $λ = 2$ with the corresponding eigenvector $[1] 0$ because

[ ] [ ] [ ] [ ] 2 1 1 = 2 = 2 1 . 0 1 0 0 0

If we rewrite the equation for an eigenvalue as

⃗ (A - λI)⃗v = 0.

We notice that this has a nonzero solution $⃗v$ only if $A - λI$ is not invertible. Were it invertible, we could write $(A - λI)-1(A - λI)⃗v = (A - λI)-1⃗0$ which implies $⃗v = ⃗0$ . Therefore, $A$ has the eigenvalue $λ$ if and only if $λ$ solves the equation

det(A - λI) = 0.

Note that this means that we will be able to find an eigenvalue without finding the corresponding eigenvector. The eigenvector will have to be found later, once $λ$ is known.

Example 3.4.2: Find all eigenvalues of $[2 1 1] 1 2 0 0 0 2$ .

We write

( ⌊ ⌋ ⌊ ⌋) ( ⌊ ⌋)
||| |||2 1 1||| |||1 0 0 |||||| ||| |||2 - λ 1 1 ||||||
det||||| |||||1 2 0|||||- λ |||||0 1 0 |||||||||| = det||||| ||||| 1 2 - λ 0 |||||||||| =
( ⌈0 0 2⌉ ⌈0 0 1 ⌉) ( ⌈ 0 0 2 - λ ⌉)
2 2
= (2 - λ )((2 - λ ) - 1) = - (λ - 1)(λ - 2)(λ - 3).

and so the eigenvalues are $λ = 1$ , $λ = 2$ , and $λ = 3$ .

Note that for an $n × n$ matrix, the polynomial we get by computing $d et(A - λI)$ will be of degree $n$ , and hence we will in general have $n$ eigenvalues.

To find an eigenvector corresponding to $λ$ , we write

(A - λI)⃗v = ⃗0,

and solve for a nontrivial (nonzero) vector $⃗v$ . If $λ$ is an eigenvalue, this will always be possible.

Example 3.4.3: Find the eigenvector of $[2 1 1] 10 2 0 02$ corresponding to the eigenvalue $λ = 3$ .

We write

(⌊ ⌋ ⌊ ⌋)⌊ ⌋ ⌊ ⌋ ⌊ ⌋ ||||2 1 1|| ||1 0 0||||||v1|| ||-1 1 1 || ||v1|| (A - λI)⃗v = ||||||||||1 2 0|||||- 3 |||||0 1 0|||||||||||||||v2|||||= |||||1 - 1 0 ||||| |||||v2|||||= ⃗0. |(|⌈ |⌉ |⌈ |⌉|)|⌈ |⌉ |⌈ |⌉ |⌈ |⌉ 0 0 2 0 0 1 v3 0 0 -1 v3

It is easy to solve this system of linear equations. Write down the augmented matrix

⌊ ⌋ ||||- 1 1 1 0|||| |||| 1 -1 0 0|||| ⌈ 0 0 - 1 0⌉

and perform row operations (exercise: which ones?) until you get

⌊|1 - 1 0 0⌋| ||||| ||||| ||⌈0 0 1 0||⌉. 0 0 0 0

The equations the entries of $⃗v$ have to satisfy are, therefore, $v1 - v2 = 0$ , $v3 = 0$ , and $v2$ is a free variable. We can pick $v2$ to be arbitrary (but nonzero) and let $v1 = v2$ and of course $v3 = 0$ . For example, $[] 1 ⃗v = 10$ . We try this:

⌊ ⌋⌊ ⌋ ⌊ ⌋ ⌊ ⌋ ||||2 1 1 ||||||1|| ||3|| ||1|| |||||1 2 0 |||||||||1|||| = ||||3|||| = 3||||1||||. ||||0 0 2 ||||||⌈ ||⌉ ||⌈ ||⌉ ||⌈ ||⌉ ⌈ ⌉ 0 0 0

Yay! It worked.

Exercise 3.4.1 (easy): Are the eigenvectors unique? Can you find a different eigenvector for $λ = 3$ in the example above? How does it relate to the other eigenvector?

Exercise 3.4.2: Note that when the matrix is $2 × 2$ you do not need to write down the augmented matrix when computing eigenvectors (if you have computed the eigenvalues correctly). Can you see why? Try it for the matrix $[2 1] 1 2$ .

3.4.2 The eigenvalue method with distinct real eigenvalues

OK. We have the equation

⃗x′ = P ⃗x.

We find the eigenvalues $λ 1$ , $λ 2$ , …, $λ n$ of the matrix $P$ , and the corresponding eigenvectors $⃗v 1$ , $⃗v 2$ , …, $⃗vn$ . Now we notice that the functions $λ1t ⃗v1e$ , $λ2t ⃗v2e$ , …, $λnt ⃗vne$ are solutions of the equation and hence $λt λ t λ t ⃗x = c1⃗v1e 1+ c2⃗v2e 2 + ⋅⋅⋅ + cn⃗vne n$ is a solution.

Theorem 3.4.1. Take $′ ⃗x = P⃗x$ . If $P$ is $n × n$ and has $n$ distinct real eigenvalues, $λ1$ , …, $λn$ then there are $n$ linearly independent corresponding eigenvectors $⃗v1$ , …, $⃗vn$ , and the general solution to the ODE can be written as.

|--------------------------------------| | ⃗x = c1⃗v1eλ1t + c2⃗v2eλ2t + ⋅⋅⋅ + cn⃗vneλnt.| ---------------------------------------

Example 3.4.4: Suppose we take the system

⌊ ⌋ |||2 1 1 ||| ⃗x′ = ||||1 2 0 ||||⃗x. |⌈0 0 2 |⌉

Find the general solution.

We have found the eigenvalues $1,2, 3$ earlier. We have found the eigenvector $[ ] 11 0$ for the eigenvalue 3. In similar fashion we find the eigenvector $[ ] -11 0$ for the eigenvalue 1 and $[ ] 01 -1$ for the eigenvalue 2 (exercise: check). Hence our general solution is

⌊ ⌋ ⌊⌋ ⌊ ⌋ ⌊ ⌋ ||| 1 ||| |||0||| |||1||| ||| et + e3t||| ⃗x = |||||- 1|||||et + |||||0||||| e2t + |||||1|||||e3t = |||||- et + e3t|||||. ⌈ 0 ⌉ ⌈1⌉ ⌈0⌉ ⌈ e2t ⌉

Exercise 3.4.3: Check that this really solves the system.

Note: If you write a homogeneous linear constant coefficient $th n$ order equation as a first order system (as we did in §3.1) then the eigenvalue equation

det(P - λI) = 0.

is essentially the same as the characteristic equation we got in §2.2 and §2.3.

3.4.3 Complex eigenvalues

A matrix might very well have complex eigenvalues even if all the entries are real. For example, suppose that we have the system

[ ] ′ 1 1 ⃗x = - 1 1 ⃗x.

Let us compute the eigenvalues of the matrix $[ ] P = 1-11 1$ .

( [ ]) 1 - λ 1 2 2 det(P - λI) = det -1 1 - λ = (1 - λ) + 1 = λ - 2λ + 2 = 0.

Thus $λ = 1 ± i$ . The corresponding eigenvectors will also be complex. First take $λ = 1 - i$ ,

⃗ ([P - (1]- i)I)⃗v = 0, i 1 -1 i⃗v = ⃗0.

It is obvious that the equations $iv1 + v2 = 0$ and $- v1 + iv2 = 0$ are multiples of each other. So we only need to consider one of them. After picking $v2 = 1$ , for example, we have the eigenvector $⃗v = [i] 1$ . In similar fashion we find that $[-i] 1$ is an eigenvector corresponding to the eigenvalue $1 + i$ .

We could write the solution as

[ ] [ ] [ (1-i)t (1+i)t] ⃗x = c1 i e(1-i)t + c2 -i e(1+i)t = c1ie(1-i)t - c2ie(1+i)t . 1 1 c1e + c2e 1

But then we would need to look for complex values $c1$ and $c2$ to solve any initial conditions. And even then it is perhaps not completely clear that we get a real solution. We could use Euler’s formula here and do the whole song and dance we did before, but we will do something a bit smarter first.

We claim that we did not have to look for the second eigenvector (nor for the second eigenvalue). All complex eigenvalues come in pairs (because the matrix $P$ is real).

First a small side note. The real part of a complex number $z$ can be computed as $z+2ˉz$ , where the bar above $z$ means $------ a + ib = a - ib$ . This operation is called the complex conjugate. Note that for a real number $a$ , $ˉa = a$ . Similarly we can bar whole vectors or matrices. If a matrix $P$ is real then $-- P = P$ . We note that $--- ---- -- P⃗x = P ⃗x = P ⃗x$ . Or

--------- -- (P - λI)⃗v = (P - λˉI)⃗v.

So if $⃗v$ is an eigenvector corresponding to eigenvalue $a + ib$ , then $-- ⃗v$ is an eigenvector corresponding to eigenvalue $a - ib$ .

Now suppose that $a + ib$ is a complex eigenvalue of $P$ , $⃗v$ the corresponding eigenvector and hence

(a+ib)t ⃗x1 = ⃗ve

is a solution (complex valued) of $⃗x′ = P ⃗x$ . Then note that $----- ea+ib = ea-ib$ and hence

-- --(a-ib)t ⃗x2 = ⃗x1 = ⃗ve

is also a solution. The function

-- (a+ib)t ⃗x1 +-⃗x1 ⃗x-1 +-⃗x2 ⃗x3 = Re ⃗x1 = R e⃗ve = 2 = 2

is also a solution. And it is real valued! Similarly as $Im z = z-ˉz- 2i$ is the imaginary part we find that

⃗x - ⃗x ⃗x4 = Im ⃗x1 =-1----2. 2i

is also a real valued solution. It turns out that $⃗x3$ and $x⃗4$ are linearly independent.

Returning to our problem, we take

[ ] [ ]( ) [ t t ] ⃗x 1 = i e(1-i)t = i etco st - iet sin t = iet cost + etsint . 1 1 e cost - ie sint

It is easy to see that

[ t ] Re ⃗x1 = e sin t , etco st [ et cost ] Im ⃗x1 = t , -e sin t

are the solutions we seek.

Exercise 3.4.4: Check that these really are solutions.

The general solution is

[et sint] [ et cost ] [c et sint + c et cos t] ⃗x = c1 t + c2 t = 1 t 2 t . e cost - e sin t c1e cos t - c2e sin t

This solution is real valued for real $c1$ and $c2$ . Now we can solve for any initial conditions that we have.

The process is this. When you have complex eigenvalues, you notice that they always come in pairs. You take one $λ = a + ib$ from the pair, you find the corresponding eigenvector $⃗v$ . You note that $(a+ib)t Re⃗ve$ and $(a+ib)t Im ⃗ve$ are also solutions to the equation, are real valued and are linearly independent. You go on to the next eigenvalue which is either a real eigenvalue or another complex eigenvalue pair. Hence, you will end up with $n$ linearly independent solutions if you had $n$ distinct eigenvalues (real or complex).

You can now find a real valued general solution to any homogeneous system where the matrix has distinct eigenvalues. When you have repeated eigenvalues, matters get a bit more complicated and we will look at that situation in §3.7.

3.4.4 Exercises

Exercise 3.4.5: Let $A$ be an $3 × 3$ matrix with an eigenvalue of 3 and a corresponding eigenvector $[ ] ⃗v = 1-1 3$ . Find $A⃗v$ .

Exercise 3.4.6: a) Find the general solution of $x′1 = 2x1$ , $x′2 = 3x2$ using the eigenvalue method (first write the system in the form $⃗x′ = A ⃗x$ ). b) Solve the system by solving each equation separately and verify you get the same general solution.

Exercise 3.4.7: Find the general solution of $x′ = 3x + x 1 1 2$ , $x′ = 2x + 4x 2 1 2$ using the eigenvalue method.

Exercise 3.4.8: Find the general solution of $x′1 = x1 - 2x2$ , $x′2 = 2x 1 + x2$ using the eigenvalue method. Do not use complex exponentials in your solution.

Exercise 3.4.9: a) Compute eigenvalues and eigenvectors of $[ ] A = -98 -23 -66 10 -2 -6$ . b) Find the general solution of $′ ⃗x = Ax⃗$ .

Exercise 3.4.10: Compute eigenvalues and eigenvectors of $[ ] -23 -12 -11 -3 -1 0$ .

3.5 Two dimensional systems and their vector fields

Note: 1 lecture, should really be in EP §5.2, but is in EP §6.2

Let us take a moment to talk about homogeneous systems in the plane. We want to think about how the vector fields look and how this depends on the eigenvalues. So we have a $2 × 2$ matrix $P$ and the system

[ ]′ [ ] x = P x . y y

(3.2)

We will be able to visually tell how the vector field looks once we find the eigenvalues and eigenvectors of the matrix.

Figure 3.3: Eigenvectors of $P$ .

Case 1. Suppose that the eigenvalues are real and positive. Find the two eigenvectors and plot them in the plane. For example, take the matrix $[1 1] 0 2$ . The eigenvalues are 1 and 2 and the corresponding eigenvectors are $[1] 0$ and $[1] 1$ . See Figure 3.3.

Now suppose that $x$ and $y$ are on the line determined by an eigenvector $⃗v$ for an eigenvalue $λ$ . That is, $[xy] = a⃗v$ for some scalar $a$ . Then

[ ]′ [ ] x = P x = P(a⃗v) = a(P ⃗v) = a λ⃗v. y y

The derivative is a multiple of $⃗v$ and hence points along the line determined by $⃗v$ . As $λ > 0$ , the derivative points in the direction of $⃗v$ when $a$ is positive and in the opposite direction when $a$ is negative. Let us draw arrows on the lines to indicate the directions. See Figure 3.4.

We fill in the rest of the arrows and we also draw a few solutions. See Figure 3.5. You will notice that the picture looks like a source with arrows coming out from the origin. Hence we call this type of picture a source or sometimes an unstable node.

Figure 3.4: Eigenvectors of $P$ with directions.

Figure 3.5: Example source vector field with eigenvectors and solutions.

Case 2. Suppose both eigenvalues were negative. For example, take the negation of the matrix in case 1, $[-1 -1] 0 -2$ . The eigenvalues are $- 1$ and $- 2$ and the corresponding eigenvectors are the same, $[1] 0$ and $[1] 1$ . The calculation and the picture are almost the same. The only difference is that the eigenvalues are negative and hence all arrows are reversed. We get the picture in Figure 3.6. We call this kind of picture a sink or sometimes a stable node.

Figure 3.6: Example sink vector field with eigenvectors and solutions.

Figure 3.7: Example saddle vector field with eigenvectors and solutions.

Case 3. Suppose one eigenvalue is positive and one is negative. For example the matrix $[1 1] 0 -2$ . The eigenvalues are 1 and $- 2$ and the corresponding eigenvectors are the same, $[1] 0$ and $[1] -3$ . We reverse the arrows on one line (corresponding to the negative eigenvalue) and we obtain the picture in Figure 3.7. We call this picture a saddle point.

The next three cases we will assume the eigenvalues are complex. In this case the eigenvectors are also complex and we cannot just plot them on the plane.

Case 4. Suppose the eigenvalues are purely imaginary. That is, suppose the eigenvalues are $± ib$ . For example, let $P = [0 1] -4 0$ . The eigenvalues turn out to be $± 2i$ and the eigenvectors are $[1] 2i$ and $[ 1] -2i$ . We take the eigenvalue $2i$ and its eigenvector $[1] 2i$ and note that the real an imaginary parts of $i2t ⃗ve$ are

[1 ] [ cos2t ] Re ei2t = , [2i] [-2 sin 2]t 1 i2t sin 2t Im 2i e = 2cos 2t .

We can take any linear combination of them, and which one we take depends on the initial conditions. For example, the real part is a a parametric equation for an ellipse. Same with the imaginary part and in fact any linear combination of them. It is not difficult to see that this is what happens in general when the eigenvalues are purely imaginary. So when the eigenvalues are purely imaginary, you get ellipses for your solutions. This type of picture is sometimes called a center. See Figure 3.8.

Figure 3.8: Example center vector field.

Figure 3.9: Example spiral source vector field.

Case 5. Now the complex eigenvalues have positive real part. That is, suppose the eigenvalues are $a ± ib$ for some $a > 0$ . For example, let $[1 1] P = -4 1$ . The eigenvalues turn out to be $1 ± 2i$ and the eigenvectors are $[1] 2i$ and $[ 1 ] -2i$ . We take $1 + 2i$ and its eigenvector $[1] 2i$ and find the real and imaginary of $⃗ve(1+2i)t$ are

[ ] [ ] Re 1 e(1+2i)t = et cos 2t , 2i -2 sin 2t [1 ] [ sin 2t ] Im e(1+2i)t = et . 2i 2cos 2t

Now note the $et$ in front of the solutions. This means that the solutions grow in magnitude while spinning around the origin. Hence we get a spiral source. See Figure 3.9.

Case 6. Finally suppose the complex eigenvalues have negative real part. That is, suppose the eigenvalues are $- a ± ib$ for some $a > 0$ . For example, let $P = [-1 -1] 4 -1$ . The eigenvalues turn out to be $- 1 ± 2i$ and the eigenvectors are $[ 1 ] -2i$ and $[1] 2i$ . We take $- 1 - 2i$ and its eigenvector $[1] 2i$ and find the real and imaginary of $(1+2i)t ⃗ve$ are

[ ] [ ] 1 (-1-2i)t -t cos2t R e 2i e = e 2 sin 2t , [ ] [ ] 1 (-1-2i)t -t- sin2t Im 2i e = e 2 cos2t .

Now note the $-t e$ in front of the solutions. This means that the solutions shrink in magnitude while spinning around the origin. Hence we get a spiral sink. See Figure 3.10.

Figure 3.10: Example spiral sink vector field.

We summarize the behavior of linear homogeneous two dimensional systems in Table 3.1.

Eigenvalues	Behavior

real and both positive	source / unstable node
real and both negative	sink / stable node
real and opposite signs	saddle
purely imaginary	center point / ellipses
complex with positive real part	spiral source
complex with negative real part	spiral sink

Table 3.1: Summary of behavior of linear homogeneous two dimensional systems.

3.5.1 Exercises

Exercise 3.5.1: Take the equation $mx ′′ + cx′ + kx = 0$ , with $m > 0$ , $c ≥ 0$ , $k > 0$ for the mass-spring system. a) Convert this to a system of first order equations. b) Classify for what $m, c,k$ do you get which behavior. c) Can you explain from physical intuition why you do not get all the different kinds of behavior here?

Exercise 3.5.2: Can you find what happens in the case when $P = [1 1] 0 1$ . In this case the eigenvalue is repeated and there is only one eigenvector. What picture does this look like?

Exercise 3.5.3: Can you find what happens in the case when $[ ] P = 11 11$ . Does this look like any of the pictures we have drawn?

3.6 Second order systems and applications

Note: more than 2 lectures, §5.3 in EP

3.6.1 Undamped mass spring systems

While we did say that we will usually only look at first order systems, it is sometimes more convenient to study the system in the way it arises naturally. For example, suppose we have 3 masses connected by springs between two walls. We could pick any higher number, and the math would be essentially the same, but for simplicity we pick 3 right now. Let us also assume no friction, that is, the system is undamped. The masses are $m1$ , $m2$ , and $m 3$ and the spring constants are $k1$ , $k2$ , $k3$ , and $k4$ . Let $x1$ be the displacement from rest position of the first mass, and $x 2$ and $x 3$ the displacement of the second and third mass. We will make, as usual, positive values go right (as $x1$ grows the first mass is moving right). See Figure 3.11.

Figure 3.11: System of masses and springs.

This simple system turns up in unexpected places. Note for example that our world really consists of small particles of matter interacting together. When we try this system with many more masses, we obtain a good approximation to how an elastic material will behave. By somehow taking a limit of the number of masses going to infinity, we obtain the continuous one dimensional wave equation. But we digress.

Let us set up the equations for the three mass system. By Hooke’s law we have that the force acting on the mass equals the spring compression times the spring constant. By Newton’s second law we have that force is mass times acceleration. So if we sum the forces acting on each mass and put the right sign in front of each term, depending on the direction in which it is acting, we end up with the desired system of equations.

m1x′1′= - k1x1 + k2(x2 - x1) = - (k1 + k2)x1 + k2x2, m x′′= - k (x - x ) + k (x - x ) = k x - (k + k )x + k x , 2 2 2 2 1 3 3 2 2 1 2 3 2 3 3 m3x′3′= - k3(x 3 - x2) - k4x3 = k3x2 - (k3 + k4)x 3.

We define the matrices

⌊|m 0 0 ⌋| ⌊|- (k + k ) k 0 ⌋| |||| 1 |||| |||| 1 2 2 |||| M = |||⌈0 m 2 0 |||⌉ and K = |||⌈ k2 - (k2 + k3) k3 |||⌉. 0 0 m3 0 k3 -(k3 + k4)

We write the equation simply as

M ⃗x ′′ = K ⃗x.

At this point we could introduce 3 new variables and write out a system of 6 equations. We claim this simple setup is easier to handle as a second order system. We will call $⃗x$ the displacement vector, $M$ the mass matrix, and $K$ the stiffness matrix.

Exercise 3.6.1: Do this setup for 4 masses (find the matrix $M$ and $K$ ). Do it for 5 masses. Can you find a prescription to do it for $n$ masses?

As before we will want to “divide by $M$ .” In this case this means computing the inverse of $M$ . All the masses are nonzero and it is easy to compute the inverse, as the matrix is diagonal.

⌊1- 0 0 ⌋ -1 |||||m1 1- ||||| M = ||||0 m2 0 ||||. ⌈0 0 1m-⌉ 3

This fact follows readily by how we multiply diagonal matrices. You should verify that $-1 -1 MM = M M = I$ as an exercise.

Let $A = M -1K$ . We look at the system $⃗x′′ = M -1K ⃗x$ , or

⃗x ′′ = A ⃗x.

Many real world systems can be modeled by this equation. For simplicity, we will only talk about the given masses-and-springs problem. We try a solution of the form

⃗x = ⃗veαt.

We compute that for this guess, $⃗x′′ = α2⃗veαt$ . We plug our guess into the equation and get

2 αt αt α ⃗ve = A⃗ve .

We can divide by $αt e$ to get that $2 α ⃗v = A⃗v$ . Hence if $2 α$ is an eigenvalue of $A$ and $⃗v$ is the corresponding eigenvector, we have found a solution.

In our example, and in many others, it turns out that $A$ has negative real eigenvalues (and possibly a zero eigenvalue). So we will study only this case here. When an eigenvalue $λ$ is negative, it means that $2 α = λ$ is negative. Hence there is some real number $ω$ such that $2 - ω = λ$ . Then $α = ±iω$ . The solution we guessed was

⃗x = ⃗v(co sωt + isin ωt).

By taking real and imaginary parts (note that $⃗v$ is real), we find that $⃗vco sωt$ and $⃗vsinωt$ are linearly independent solutions.

If an eigenvalue is zero, it turns out that $⃗v$ and $⃗vt$ are solutions if $⃗v$ is the corresponding eigenvector.

Exercise 3.6.2: Show that if $A$ has a zero eigenvalue and $⃗v$ is the corresponding eigenvector, then $⃗x =⃗v(a + bt)$ is a solution of $⃗x ′′ = A ⃗x$ for arbitrary constants $a$ and $b$ .

Theorem 3.6.1. Let $A$ be an $n × n$ with $n$ distinct real negative eigenvalues we denote by $2 2 2 - ω1 > -ω 2 > ⋅⋅⋅ > -ω n$ , and the corresponding eigenvectors by $⃗v1$ , $⃗v2$ , …, $⃗vn$ . If $A$ is invertible (that is, if $ω1 > 0$ ), then

|----------------------------------| | ∑n | | ⃗x(t) = ⃗vi(aicos ωit + bisinωit),| --------i=1--------------------------

is the general solution of

′′ ⃗x = A ⃗x,

for some arbitrary constants $ai$ and $bi$ . If $A$ has a zero eigenvalue, that is $ω 1 = 0$ , and all other eigenvalues are distinct and negative then the general solution becomes

|--------------------∑n--------------------------| | ⃗x(t) = ⃗v (a + b t) + ⃗v(a cosω t + b sin ω t). | | 1 1 1 i i i i i | ----------------------i=2-------------------------|

Note that we can use this solution and the setup from the introduction of this section even when some of the masses and springs are missing. For example, when there are say 2 masses and only 2 springs, simply take only the equations for the two masses and set all the spring constants for the springs that are missing to zero.

3.6.2 Examples

Example 3.6.1: Suppose we have the system in Figure 3.12, with $m 1 = 2$ , $m 2 = 1$ , $k1 = 4$ , and $k2 = 2$ .

Figure 3.12: System of masses and springs.

The equations we write down are

[ ] [ ] 2 0 ⃗x′′ = - (4 + 2) 2 ⃗x, 0 1 2 -2

[ ] ′′ -3 1 ⃗x = 2 - 2 ⃗x.

We find the eigenvalues of $A$ to be $λ = - 1,-4$ (exercise). Now we find the eigenvectors to be $[1] 2$ and $[1] -1$ respectively (exercise).

We check the theorem and note that $ω 1 = 1$ and $ω 2 = 2$ . Hence the general solution is

[ ] [ ] ⃗x = 1 (a co st + b sint) + 1 (a cos2t + b sin 2t). 2 1 1 -1 2 2

The two terms in the solution represent the two so-called natural or normal modes of oscillation. And the two (angular) frequencies are the natural frequencies. The two modes are plotted in Figure 3.13.

Figure 3.13: The two modes of the mass spring system. In the left plot the masses are moving in unison and the right plot are masses moving in the opposite direction.

Let us write the solution as

[ ] [ ] 1 1 ⃗x = 2 c 1cos(t - α2) + - 1 c2co s(2t - α1).

The first term,

[1] [ c co s(t - α )] c1cos(t - α1) = 1 1 , 2 2c1cos(t - α1)

corresponds to the mode where the masses move synchronously in the same direction.

On the other hand the second term,

[ ] [ ] 1 c2cos(2t - α 2) - 1 c2cos(2t - α 2) = - c2cos(2t - α2) ,

corresponds to the mode where the masses move synchronously but in opposite directions.

The general solution is a combination of the two modes. That is, the initial conditions determine the amplitude and phase shift of each mode.

Example 3.6.2: Let us do another example. In this example we have two toy rail cars. Car 1 of mass 2 kg is traveling at 3 ^m/_s towards the second rail car of mass 1 kg. There is a bumper on the second rail car which engages one the cars hit (it connects to two cars) and does not let go. The bumper acts like a spring of spring constant $k = 2N∕m$ . The second car is 10 meters from a wall. See Figure 3.14.

Figure 3.14: The crash of two rail cars.

We want to ask several question. At what time after the cars link does impact with the wall happen? What is the speed of car 2 when it hits the wall?

OK, let us first set the system up. Let us assume that time $t = 0$ is the time when the two cars link up. Let $x1$ be the displacement of the first car from the position at $t = 0$ , and let $x2$ be the displacement of the second car from its original location. Then the time when $x (t) = 10 2$ is exactly the time when impact with wall occurs. For this $t$ , $′ x2(t)$ is the speed at impact. This system acts just like the system of the previous example but without $k1$ . Hence the equation is

[ ] [ ] 2 0 ⃗x′′ = -2 2 ⃗x. 0 1 2 - 2

[ ] ⃗x ′′ = -1 1 ⃗x. 2 - 2

We compute the eigenvalues of $A$ . It is not hard to see that the eigenvalues are 0 and $- 3$ (exercise). Furthermore, the eigenvectors are $[1] 1$ and $[1] -2$ respectively (exercise). We note that $√-- ω 2 = 3$ and we use the second part of the theorem to find our general solution to be

[1] [ 1 ]( √ -- √ -)
⃗x = (a1 + b 1t) + a2 cos 3t + b2 sin 3t =
1 - 2 -- --
[ a1 + b1t + a2cos √ 3t + b2sin √3t ]
= √ -- √--
a1 + b 1t - 2a2cos 3t - 2b2sin 3t

We now apply the initial conditions. First the cars start at position 0 so $x1(0) = 0$ and $x2(0) = 0$ . The first car is traveling at 3 ^m/_s, so $x′1(0) = 3$ and the second car starts at rest, so $x′2(0 ) = 0$ . The first conditions says

[ ] ⃗ a1 + a2 0 = ⃗x(0) = a1 - 2a2 .

It is not hard to see that this implies that $a1 = a 2 = 0$ . We plug $a1$ and $a 2$ and differentiate to get

[ √-- √-- ] ′ b1 + 3b 2cos 3t ⃗x(t) = √ -- √ -- . b 1 - 2 3b2 cos 3t

[3 ] [ b + √3b ] = ⃗x′(0) = 1 √ --2 . 0 b 1 - 2 3b2

It is not hard to solve these two equations to find $b1 = 2$ and $b 2 =-1√- 3$ . Hence the position of our cars is (until the impact with the wall)

⌊ 1 √--⌋ ||||2t + √3 sin √3t|||| ⃗x = |⌈2t - 2√-sin 3t|⌉. 3

Note how the presence of the zero eigenvalue resulted in a term containing $t$ . This means that the carts will be traveling in the positive direction as time grows, which is what we expect.

What we are really interested in is the second expression, the one for $x2$ . We have $2√-- √-- x2(t) = 2t - 3 sin 3t$ . See Figure 3.15 for the plot of $x2$ versus time.

Figure 3.15: Position of the second car in time (ignoring the wall).

Just from the graph we can see that time of impact will be a little more than 5 seconds from time zero. For this you have to solve the equation $2-- √-- 10 = x2(t) = 2t - √3 sin 3t$ . Using a computer (or even a graphing calculator) we find that $timpact ≈ 5.22$ seconds.

As for the speed we note that $√ -- x′2 = 2 - 2 cos 3t$ . At time of impact (5.22 seconds from $t = 0$ ) we get that $x′(timpact) ≈ 3.8 5 2$ .

The maximum speed is the maximum of $√ -- 2 - 2cos 3t$ , which is 4. We are traveling at almost the maximum speed when we hit the wall.

Now suppose that Bob is a tiny person sitting on car 2. Bob has a Martini in his hand and would like to not spill it. Let us suppose Bob would not spill his Martini when the first car links up with car 2, but if car 2 hits the wall at any speed greater than zero, Bob will spill his drink. Suppose Bob can move car 2 a few meters towards or away from the wall (he cannot go all the way to the wall, nor can he get out of the way of the first car). Is there a “safe” distance for him to be in? A distance such that the impact with the wall is at zero speed?

The answer is yes. Looking at Figure 3.15, we note the “plateau” between $t = 3$ and $t = 4$ . There is a point where the speed is zero. To find it we need to solve $x′2(t) = 0$ . This is when $-- cos √ 3t = 1$ or in other words when $t = 2√π, 4√π,... 3 3$ and so on. We plug in the first value to obtain $( ) x2 2√π3- = 4√π3-≈ 7.26$ . So a “safe” distance is about 7 and a quarter meters from the wall.

Alternatively Bob could move away from the wall towards the incoming car 2 where another safe distance is $8π- √3 ≈ 14.51$ and so on, using all the different $t$ such that $′ x2(t) = 0$ . Of course $t = 0$ is always a solution here, corresponding to $x2 = 0$ , but that means standing right at the wall.

3.6.3 Forced oscillations

Finally we move to forced oscillations. Suppose that now our system is

′′ ⃗x = A ⃗x +F⃗co sωt.

(3.3)

That is, we are adding periodic forcing to the system in the direction of the vector $F⃗$ .

Just like before this system just requires us to find one particular solution $⃗xp$ , add it to the general solution of the associated homogeneous system $⃗x c$ and we will have the general solution to (3.3). Let us suppose that $ω$ is not one of the natural frequencies of $′′ ⃗x = A⃗x$ , then we can guess

⃗x = ⃗cco sωt, p

where $⃗c$ is an unknown constant vector. Note that we do not need to use sine since there are only second derivatives. We solve for $⃗c$ to find $⃗xp$ . This is really just the method of undetermined coefficients for systems. Let us differentiate $⃗xp$ twice to get

′′ 2 ⃗xp = -ω ⃗ccosωt.

Now plug into the equation

- ω2⃗c cosωt = A⃗ccosωt + ⃗F cosωt

We can cancel the cosine and rearrange to obtain

2 ⃗ (A + ω I)⃗c = -F.

⃗c = (A + ω 2I)-1(-F⃗).

Of course this means that $2 2 (A + ω I) = (A - (-ω )I)$ is invertible. That matrix is invertible if and only if $2 - ω$ is not an eigenvalue of $A$ . That is true if and only if $ω$ is not a natural frequency of the system.

Example 3.6.3: Let us take the example in Figure 3.12 with the same parameters as before: $m 1 = 2$ , $m 2 = 1$ , $k1 = 4$ , and $k2 = 2$ . Now suppose that there is a force $2 cos3t$ acting on the second cart.

The equation is

[ ] [ ] ′′ -3 1 0 ⃗x = 2 -2 ⃗x + 2 cos3t.

We have solved the associated homogeneous equation before and found the complementary solution to be

[] [ ] 1 1 ⃗xc = 2 (a 1cost + b 1sin t) + -1 (a2cos2t + b2sin2t).

We note that the natural frequencies were 1 and 2. Hence 3 is not a natural frequency, we can try $⃗c cos3t$ . We can invert $(A + 32I)$

([- 3 1 ] )-1 [6 1]-1 [-7 -1] + 32I = = 4-01 403 . 2 -2 2 7 20 20

Hence,

[ ] [ ] [ ] 2 -1 470 -410 0 210 ⃗c = (A + ω I) (- ⃗F) = -1 -3 = -3 . 20 20 -2 10

Combining with what we know the general solution of the associated homogeneous problem to be we get that the general solution to $⃗x′′ = A⃗x + ⃗F cosωt$ is

[1] [ 1 ] [ 1-] ⃗x = ⃗xc + ⃗xp = (a1cos t + b1sint) + (a2cos 2t + b2sin2t) + 2-03 cos3t. 2 - 1 10

The constants $a 1$ , $a2$ , $b1$ , and $b2$ must then be solved for given any initial conditions.

If $ω$ is a natural frequency of the system resonance occurs because you will have to try a particular solution of the form

⃗xp = ⃗ctsinωt + ⃗d cosωt.

That is assuming that all eigenvalues of the coefficient matrix are distinct. Note that the amplitude of this solution grows without bound as $t$ grows.

3.6.4 Exercises

Exercise 3.6.3: Find a particular solution to

[ ] [ ] ⃗x′′ = -3 1 ⃗x + 0 cos2t. 2 -2 2

Exercise 3.6.4: Let us take the example in Figure 3.12 with the same parameters as before: $m 1 = 2$ , $k1 = 4$ , and $k2 = 2$ , except for $m 2$ which is unknown. Suppose that there is a force $cos5t$ acting on the first mass. Find an $m 1$ such that there exists a particular solution where the first mass does not move.

Note: This idea is called dynamic damping. In practice there will be a small amount of damping and so any transient solution will disappear and after long enough time, the first mass will always come to a stop.

Exercise 3.6.5: Let us take the example 3.6.2, but that at time of impact, cart 2 is moving to the left at the speed of 3 ^m/_s. a) Find the behavior of the system after linkup. b) Will the second car hit the wall, or will it be moving away from the wall as time goes on. c) at what speed would the first car have to be traveling for the system to essentially stay in place after linkup.

Exercise 3.6.6: Let us take the example in Figure 3.12 with parameters $m 1 = m2 = 1$ , $k1 = k2 = 1$ . Does there exist a set of initial conditions for which the first cart moves but the second cart does not? If so find those conditions, if not argue why not.

3.7 Multiple eigenvalues

Note: 1–2 lectures, §5.4 in EP

It may very well happen that a matrix has some “repeated” eigenvalues. That is, the characteristic equation $det(A - λI) = 0$ may have repeated roots. As we have said before, this is actually unlikely to happen for a random matrix. If you take a small perturbation of $A$ (you change the entries of $A$ slightly) you will get a matrix with distinct eigenvalues. As any system you will want to solve in practice is an approximation to reality anyway, it is not indispensable to know how to solve these corner cases. But it may happen on occasion that it is easier or desirable to solve such a system directly.

3.7.1 Geometric multiplicity

Take the diagonal matrix

[ ] 3 0 A = 0 3 .

$A$ has an eigenvalue 3 of multiplicity 2. We usually call the multiplicity of the eigenvalue in the characteristic equation the algebraic multiplicity. In this case, there exist 2 linearly independent eigenvectors, $[1] 0$ and $[0] 1$ . This means that the so-called geometric multiplicity of this eigenvalue is 2.

In all the theorems where we required a matrix to have $n$ distinct eigenvalues, we only really needed to have $n$ linearly independent eigenvectors. For example, let $⃗x′ = A ⃗x$ has the general solution

[ ] [ ] ⃗x = c 1 e3t + c 0 e3t. 1 0 2 1

Let us restate the theorem about real eigenvalues. In the following theorem we will repeat eigenvalues according to (algebraic) multiplicity. So for $A$ above we would say that it has eigenvalues 3 and 3.

Theorem 3.7.1. Take $⃗x′ = P⃗x$ . If $P$ is $n × n$ and has $n$ real eigenvalues (not necessarily distinct), $λ1$ , …, $λn$ , and if there are $n$ linearly independent corresponding eigenvectors $⃗v1$ , …, $⃗vn$ , and the general solution to the ODE can be written as.

⃗x = c1⃗v1eλ1t + c2⃗v2eλ2t + ⋅⋅⋅ + cn⃗vneλnt.

The geometric multiplicity of an eigenvalue of algebraic multiplicity $n$ is equal to the number of linearly independent eigenvectors we can find. It is not hard to see that the geometric multiplicity is always less than or equal to the algebraic multiplicity. Above we, therefore, handled the case when these two numbers are equal. If the geometric multiplicity is equal to the algebraic multiplicity we say the eigenvalue is complete.

The hypothesis of the theorem could, therefore, be stated as saying that if all the eigenvalues of $P$ are complete then there are $n$ linearly independent eigenvectors and thus we have the given general solution.

Note that if the geometric multiplicity of an eigenvalue is 2 or greater, then the set of linearly independent eigenvectors is not unique up to multiples as it was before. For example, for the diagonal matrix $A$ above we could also pick eigenvectors $[1] 1$ and $[1 ] -1$ , or in fact any pair of two linearly independent vectors.

3.7.2 Defective eigenvalues

If an $n × n$ matrix has less than $n$ linearly independent eigenvectors, it is said to be deficient. Then there is at least one eigenvalue with algebraic multiplicity that is higher than the geometric multiplicity. We call this eigenvalue defective and the difference between the two multiplicities we call the defect.

Example 3.7.1: The matrix

[ ] 3 1 0 3

has an eigenvalue 3 of algebraic multiplicity 2. Let us try to compute the eigenvectors.

[0 1] [v ] 1 = ⃗0. 0 0 v2

We must have that $v2 = 0$ . Hence any eigenvector is of the form $[v01]$ . Any two such vectors are linearly dependent, and hence the geometric multiplicity of the eigenvalue is 1. Therefore, the defect is 1, and we can no longer apply the eigenvalue method directly to a system of ODEs with such a coefficient matrix.

The key observation we will use here is that if $λ$ is an eigenvalue of $A$ of algebraic multiplicity $m$ , then we will be able to find $m$ linearly independent vectors solving the equation $(A - λI)m⃗v = ⃗0$ . We will call these the generalized eigenvectors.

Let us continue with the example $[3 1] A = 0 3$ and the equation $′ ⃗x = A ⃗x$ . We have an eigenvalue $λ = 3$ of (algebraic) multiplicity 2 and defect 1. We have found one eigenvector $[ ] ⃗v1 = 10$ . We have the solution

3t ⃗x1 = ⃗v1e .

In this case, let us try (in the spirit of repeated roots of the characteristic equation for a single equation) another solution of the form

3t ⃗x2 = (⃗v2 + ⃗v1t)e .

We differentiate to get

′ 3t 3t 3t 3t ⃗x2 = ⃗v1e + 3(⃗v2 + ⃗v1t)e = (3⃗v2 + ⃗v1)e + 3⃗v1te .

$⃗x2′$ must equal $A ⃗x2$ , and

A⃗x2 = A(⃗v2 + ⃗v1t)e3t = A⃗v2e3t + A⃗v1te3t.

By looking at the coefficients of $e3t$ and $te3t$ we see $3⃗v2 + ⃗v1 = A ⃗v2$ and $3⃗v1 = A⃗v1$ . This means that

⃗ (A - 3I)⃗v1 = 0, and (A - 3I)⃗v2 = ⃗v1.

If these two equations are satisfied, then $x⃗2$ is a solution. We know the first of these equations is satisfied because $⃗v1$ is an eigenvector. If we plug the second equation into the first we find that

(A - 3I)(A - 3I)⃗v2 = ⃗0, or (A - 3I)2⃗v2 = ⃗0.

If we can, therefore, find a $⃗v2$ which solves $(A - 3I)2⃗v2 = ⃗0$ , and such that $(A - 3I)⃗v2 = ⃗v1$ , we are done. This is just a bunch of linear equations to solve and we are by now very good at that.

We notice that in this simple case $(A - 3I)2$ is just the zero matrix (exercise). Hence, any vector $⃗v2$ solves $(A - 3I)2⃗v2 = ⃗0$ . So we just have to make sure that $(A - 3I)⃗v2 = ⃗v1$ . Write

[ ] [ ] [ ] 0 1 a = 1 . 0 0 b 0

By inspection we see that letting $a = 0$ ( $a$ could be anything in fact) and $b = 1$ does the job. Hence we can take $⃗v = [0] 2 1$ . So our general solution to $⃗x′ = A ⃗x$ is

[ ] ([ ] [ ] ) [ ] 1 3t 0 1 3t c1e3t + c2te3t ⃗x = c1 0 e + c2 1 + 0 t e = c e3t . 2

Let us check that we really do have the solution. First $′ 3t 3t 3t x 1 = c13e + c2e + 3c2te = 3x 1 + x2$ , good. Now $x′2 = 3c2e3t = 3x2$ , good.

Note that the system $′ ⃗x = A ⃗x$ has a simpler solution since $A$ is a triangular matrix. In particular, the equation for $x2$ does not depend on $x1$ .

Exercise 3.7.1: Solve $[ ] ⃗x′ = 30 13 ⃗x$ by first solving for $x2$ and then for $x 1$ independently. Now check that you got the same solution as we did above.

Let us describe the general algorithm. First for $λ$ of multiplicity 2, defect 1. First find an eigenvector $⃗v1$ of $λ$ . Now find a vector $⃗v2$ such that

2 (A - 3I)⃗v2 = ⃗0, (A - 3I)⃗v2 = ⃗v1.

This gives us two linearly independent solutions

λt ⃗x1 = ⃗v(1e , ) ⃗x2 = ⃗v2 + ⃗v1t eλt.

This machinery can also be generalized to larger matrices and higher defects. We will not go over, but let us just state the ideas. Suppose that $A$ has a multiplicity $m$ eigenvalue $λ$ . We find vectors such that

(A - λI)k⃗v = ⃗0, but (A - λI)k-1⃗v ⇔ ⃗0.

Such vectors are called generalized eigenvectors. For every eigenvector $⃗v1$ we find a chain of generalized eigenvectors $⃗v2$ through $⃗vk$ such that:

(A - λI)⃗v1 = ⃗0, (A - λI)⃗v2 =⃗v1, ... (A - λI)⃗vk =⃗vk-1.

We form the linearly independent solutions

λt ⃗x1 = ⃗v1e , ⃗x2 = (⃗v2 + ⃗v1t)eλt, ... ( k-2 k-1 ) ⃗x = ⃗v + ⃗v t + ⋅⋅⋅ + ⃗v-t----+ ⃗v -t----- eλt. k k k-1 2(k - 2)! 1 (k - 1)!

We proceed to find chains until we form $m$ linearly independent solutions ( $m$ is the multiplicity). You may need to find several chains for every eigenvalue.

3.7.3 Exercises

Exercise 3.7.2: Let $[ ] A = 53 - -31$ . Solve $′ ⃗x = A⃗x$ .

Exercise 3.7.3: Let $[ 5 -4 4 ] A = 0 3 0 -2 4 -1$ . a) What are the eigenvalues? b) What is/are the defect(s) of the eigenvalue(s)? c) Solve $′ ⃗x = A⃗x$ .

Exercise 3.7.4: Let $[ ] A = 20 1 2 00 0 0 2$ . a) What are the eigenvalues? b) What is/are the defect(s) of the eigenvalue(s)? c) Solve $′ ⃗x = A⃗x$ in two different ways and verify you get the same answer.

Exercise 3.7.5: Let $[ ] -01 -12 2-2 A = -4 4 7$ . a) What are the eigenvalues? b) What is/are the defect(s) of the eigenvalue(s)? c) Solve $⃗x ′ = A⃗x$ .

Exercise 3.7.6: Let $[ 0 4 -2] A = -10 -40 1-2$ . a) What are the eigenvalues? b) What is/are the defect(s) of the eigenvalue(s)? c) Solve $⃗x ′ = A⃗x$ .

Exercise 3.7.7: Let $[ 2 1 -1] A = -1 0 2 -1 -2 4$ . a) What are the eigenvalues? b) What is/are the defect(s) of the eigenvalue(s)? c) Solve $′ ⃗x = A⃗x$ .

Exercise 3.7.8: Suppose that $A$ is a $2 × 2$ matrix with a repeated eigenvalue $λ$ . Suppose that there are two linearly independent eigenvectors. Show that the matrix is diagonal, in particular $A = λI$ .

3.8 Matrix exponentials

Note: 2 lectures, §5.5 in EP

3.8.1 Definition

In this section we present a different way of finding the fundamental matrix solution of a system. Suppose that we have the constant coefficient equation

⃗x′ = Px⃗,

as usual. Now suppose that this was one equation ( $P$ is a number or a $1 × 1$ matrix). Then the solution to this would be

⃗x = ePt.

It turns out the same computation works for matrices when we define $ePt$ properly. First let us write down the Taylor series for $eat$ for some number $a$ .

2 3 4 ∑∞ k at (at)-- (at)-- (at)- (at)-- e = 1 + at + 2 + 6 + 24 + ⋅⋅⋅ = k! . k=0

Recall $k! = 1 ⋅ 2 ⋅ 3⋅⋅⋅k$ , and $0! = 1$ . We differentiate this series term by term

32 4 3 ( 2 3 ) a + a2t + a-t-+ a-t-+ ⋅⋅⋅ = a 1 + at + (at)-+ (at)-+ ⋅⋅⋅ = aeat. 2 6 2 6

Maybe we can write try the same trick here. Suppose that for an $n × n$ matrix $A$ we define the matrix exponential as

|------------------------------------------| |eA d=efI + A + 1-A2 + 1A 3 + ⋅⋅⋅ + 1-Ak + ⋅⋅⋅ --------------2------6----------k!----------

Let us not worry about convergence. The series really does always converge. We usually write $Pt$ as $tP$ by convention when $P$ is a matrix. With this small change and by the exact same calculation as above we have that

d ( tP ) tP -- e = Pe . dt

Now $P$ and hence $tP e$ is an $n × n$ matrix. What we are looking for is a vector. We note that in the $1 × 1$ case we would at this point multiply by an arbitrary constant to get the general solution. In the matrix case we multiply by a column vector $⃗c$ .

Theorem 3.8.1. Let $P$ be an $n × n$ matrix. Then the general solution to $′ ⃗x = P ⃗x$ is

tP ⃗x = e ⃗c,

where $⃗c$ is an arbitrary constant vector. In fact $⃗x(0) = ⃗c$ .

Let us check.

( ) d-⃗x = d- etP⃗c = PetP⃗c = P ⃗x. dt dt

Hence $etP$ is the fundamental matrix solution of the homogeneous system. If we find a way to compute the matrix exponential, we will have another method of solving constant coefficient homogeneous systems. It also makes it easy to solve for initial conditions. To solve $′ ⃗x = Ax⃗$ , $⃗ ⃗x(0) = b$ , we take the solution

⃗x = etA⃗b.

This equation follows because $e0A = I$ , so $⃗x(0) = e0A⃗b = ⃗b$ .

We mention a drawback of matrix exponentials. In general $A+B A B e ⇔ e e$ . The trouble is that matrices do not commute, that is, in general $AB ⇔ BA$ . If you try to prove $eA+B ⇔ eAeB$ using the Taylor series, you will see why the lack of commutativity becomes a problem. However, it is still true that if $AB = BA$ , that is, if $A$ and $B$ commute, then $eA+B = eAeB$ . We will find this fact useful. Let us restate this as a theorem to make a point.

Theorem 3.8.2. If $AB = BA$ then $A+B A B e = e e$ . Otherwise $A+B A B e ⇔ e e$ in general.

3.8.2 Simple cases

In some instances it may work to just plug into the series definition. Suppose the matrix is diagonal. For example, $D = [a 0] 0 b$ . Then

[ k ] Dk = a 0k , 0 b

and

[ ] [ ] [ 2 ] [ 3 ] [ a ] eD = I + D + 1-D2 + 1D 3 + ⋅⋅⋅ = 1 0 + a 0 + 1- a 02 + 1- a 03 + ⋅⋅⋅ = e 0b . 2 6 0 1 0 b 2 0 b 6 0 b 0 e

So by this rationale we have that

[ ] [ a ] eI = e 0 and eaI = e 0 . 0 e 0 ea

This makes exponentials of certain other matrices easy to compute. Notice for example that the matrix $[ ] A = 53 --31$ can be written as $2I + B$ where $[ ] B = 33 - -33$ . Notice that $2I$ and $B$ commute, and that $B 2 = [00 0 0]$ . So $Bk = 0$ for all $k ≥ 2$ . Therefore, $eB = I + B$ . Suppose we actually want to compute $etA$ . $2tI$ and $tB$ still commute (exercise: check this) and $etB = I + tB$ , since $(tB)2 = t2B 2 = 0$ . We write

[ ]
tA 2tI+tB 2tI tB e2t 0
e = e = e e = 0 e2t (I + tB ) =
[ 2t ] [ ] [ 2t 2t ]
= e 0 1 + 3t -3t = (1 + 3t)e - 3te .
0 e2t 3t 1 - 3t 3te2t (1 - 3t)e2t

So we have found the fundamental matrix solution for the system $′ ⃗x = A⃗x$ . Note that this matrix has a repeated eigenvalue with a defect; there is only one eigenvector for the eigenvalue 2. So we have found a perhaps easier way to handle this case. In fact, if a matrix $A$ is $2 × 2$ and has an eigenvalue $λ$ of multiplicity 2, then either it is diagonal, or $A = λI + B$ where $B2 = 0$ . This is a good exercise.

Exercise 3.8.1: Suppose that $A$ is $2 × 2$ and $λ$ is the only eigenvalue. Then show that $2 (A - λI) = 0$ . Then we can write $A = λI + B$ , where $2 B = 0$ . Hint: First write down what does it mean for the eigenvalue to be of multiplicity 2. You will get an equation for the entries. Now compute the square of $B$ .

Matrices $B$ such that $Bk = 0$ for some $k$ are called nilpotent. Computation of the matrix exponential for nilpotent matrices is easy by just writing down the first $k$ terms of the Taylor series.

3.8.3 General matrices

In general, the exponential is not as easy to compute as above. We cannot usually write any matrix as a sum of commuting matrices where the exponential is simple for each one. But fear not, it is still not too difficult provided we can find enough eigenvectors. First we need the following interesting result about matrix exponentials. For any two square matrices $A$ and $B$ , we have

eBAB-1 = BeAB -1.

This can be seen by writing down the Taylor series. First note that

-1 2 -1 -1 -1 2 -1 (BAB ) = BAB BAB = BAIAB = BA B .

And hence by the same reasoning $-1k k -1 (BAB ) = BA B$ . So now write down the Taylor series for $-1 eBAB$

eBAB-1 = I + BAB -1 + 1(BAB -1)2 + 1-(BAB -1)3 + ⋅⋅⋅ 2 6 -1 -1 1- 2 -1 1- 3 -1 = BB + BAB + 2 BA B + 6BA B + ⋅⋅⋅ ( 1 1 ) = B I + A +--A2 + -A3 + ⋅⋅⋅B -1 2 6 = BeAB -1.

Now we will write a general matrix $A$ as $EDE -1$ , where $D$ is diagonal. This procedure is called diagonalization. If we can do that, you can see that the computation of the exponential becomes easy. Adding $t$ into the mix we see that

tA tD -1 e = Ee E .

Now to do this we will need $n$ linearly independent eigenvectors of $A$ . Otherwise this method does not work and we need to be trickier, but we will not get into such details in this course. We let $E$ be the matrix with the eigenvectors as columns. Let $λ1$ , …, $λn$ be the eigenvalues and let $⃗v1$ , …, $⃗vn$ be the eigenvectors, then $E = [⃗v1 ⃗v2 ⋅⋅⋅ ⃗vn]$ . Let $D$ be the diagonal matrix with the eigenvalues on the main diagonal. That is

⌊ ⌋ ||||λ1 0 ⋅⋅⋅ 0 |||| |||| 0 λ2 ⋅⋅⋅ 0 |||| D = |||| .. .. .. ..||||. |||⌈ . . . .|||⌉ 0 0 ⋅⋅⋅ λn

Now we write

AE = A[⃗v1 ⃗v2 ⋅⋅⋅ ⃗vn] = [A ⃗v1 A⃗v2 ⋅⋅⋅ A ⃗vn] = [λ 1⃗v1 λ2⃗v2 ⋅⋅⋅ λn⃗vn] = [⃗v ⃗v ⋅⋅⋅ ⃗v ]D 1 2 n = ED.

Now the columns of $E$ are linearly independent as these are the eigenvectors of $A$ . Hence $E$ is invertible. Since $AE = ED$ , we right multiply by $E -1$ and we get

A = EDE -1.

This means that $eA = EeDE -1$ . With $t$ is turns into

|--------------------------------------------| | ⌊λ1t ⌋ | | |||||e 0λ2t ⋅⋅⋅ 0 ||||| | | tA tD -1 ||||0 e ⋅⋅⋅ 0 |||| -1 | | e = Ee E = E |||| ... ... ... ... ||||E . | | ||⌈ λnt||⌉ | ---------------------0----0---⋅⋅⋅--e---------

(3.4)

The formula (3.4), therefore, gives the formula for computing the fundamental matrix solution $etA$ for the system $′ ⃗x = A⃗x$ , in the case where we have $n$ linearly independent eigenvectors.

Notice that this computation still works when the eigenvalues and eigenvectors are complex, though then you will have to compute with complex numbers. Note that it is clear from the definition that if $A$ is real, then $etA$ is real. So you will only need complex numbers in the computation and you may need to apply Euler’s formula to simplify the result. If simplified properly the final matrix will not have any complex numbers in it.

Example 3.8.1: Compute the fundamental matrix solution using the matrix exponentials for the system

[ ]′ [ ] [ ] x = 1 2 x . y 2 1 y

Then compute the particular solution for the initial conditions $x(0) = 4$ and $y(0) = 2$ .

Let $A$ be the coefficient matrix $[1 2] 2 1$ . We first compute (exercise) that the eigenvalues are 3 and $- 1$ and the corresponding eigenvectors are $[1] 1$ and $[1 ] -1$ . Hence we write

[1 1 ] [e3t 0 ] [1 1 ]-1 etA = -t 1[ - 1] [0 e ] 1 [-1 ] 1 1 e3t 0 -1- -1 - 1 = 1 - 1 0 e-t 2 -1 1 [3t -t ] [ ] = --1 e e - 1 -1 2 e3t -e-t - 1 1 [ 3t -t 3t -t] [e3t+e-t e3t-e-t] = --1 - e - e - e + e = 3t2-t 3t2-t . 2 - e3t + e-t - e3t - e-t e--2e- e+e2--

The initial conditions are $x(0) = 4$ and $y(0) = 2$ . Hence, by the property that $e0A = I$ we find that the particular solution we are looking for is $etA⃗b$ where $⃗b$ is $[4] 2$ . Then the particular solution we are looking for is

[ ] [ 3t -t 3t -t] [ ] [ ] [ ] x e-+2e- e--2e- 4 2e3t + 2e-t + e3t - e-t 3e3t + e-t y = e3t-e-t e3t+e-t 2 = 2e3t - 2e-t + e3t + e-t = 3e3t - e-t . 2 2

3.8.4 Fundamental matrix solutions

We note that if you can compute the fundamental matrix solution in a different way, you can use this to find the matrix exponential $tA e$ . The fundamental matrix solution of a system of ODEs is not unique. The exponential is the fundamental matrix solution with the property that for $t = 0$ we get the identity matrix. So we must find the right fundamental matrix solution. Let $X$ be any fundamental matrix solution to $⃗x′ = A⃗x$ . Then we claim

tA -1 e = X(t)[X(0)] .

Obviously if we plug $t = 0$ into $X(t)[X (0)]-1$ we get the identity. It is not hard to see that we can multiply a fundamental matrix solution on the right by any constant invertible matrix and we still get a fundamental matrix solution. All we are doing is changing what the arbitrary constants are in the general solution $⃗x(t) = X (t)⃗c$ .

3.8.5 Approximations

If you think about it, the computation of any fundamental matrix solution $X$ using the eigenvalue method is just as difficult as computation of $etA$ . So perhaps we did not gain much by this new tool. However, the Taylor series expansion actually gives us a very easy way to approximate solutions, which the eigenvalue method did not.

The simplest thing we can do is to just compute the series up to a certain number of terms. There are better ways to approximate the exponential . In many cases however, few terms of the Taylor series give a reasonable approximation for the exponential and may suffice for the application. For example, let us compute the first 4 terms of the series for the matrix $[1 2] A = 2 1$ .

[ ] [ ] [ ]
tA t2 2 t3 3 1 2 2 52 2 3 163 73
e ≈ I + tA + 2 A + 6 A = I + t 2 1 + t 2 5 + t 7 13 =
2 [ 3 6 ]
1 + t + 52t2 + 163t3 2t + 2t2 + 73t3
= 2 7 3 5 2 13 3 .
2t + 2t + 3t 1 + t + 2t + 6 t

Just like the Taylor series approximation for the scalar version, the approximation will be better for small $t$ and worse for larger $t$ . For larger $t$ , you will generally have to compute more terms. Let us see how we stack up against the real solution with $t = 0.1$ . The approximate solution is approximately (rounded to 8 decimal places)

⌊ ⌋ 2 3 ||1.127166 67 0.2 2233333 || e0.1A ≈ I + 0.1A + 0.1-A2 + 0.1-A3 = |||||0.222333 33 1.1 2716667 |||||. 2 6 |⌈ |⌉

And plugging $t = 0.1$ into the real solution (rounded to 8 decimal places) we get

[ ] e0.1A = 1.1273481 1 0.22 251069 . 0.2225106 9 1.12 734811

This is not bad at all. Although if you take the same approximation for $t = 1$ you get (using the Taylor series)

[6.6666 6667 6.333333 33] , 6.3333 3333 6.666666 67

while the real value is (again rounded to 8 decimal places)

[ ] 10.226 70818 9.8588287 4 9.8588287 4 10.226 70818 .

So the approximation is not very good once we get up to $t = 1$ . To get a good approximation at $t = 1$ (say up to 2 decimal places) you would need to go up to the $11th$ power (exercise).

3.8.6 Exercises

Exercise 3.8.2: Find a fundamental matrix solution for the system $x′ = 3x + y$ , $y′ = x + 3y$ .

Exercise 3.8.3: Find $eAt$ for the matrix $[ ] A = 20 3 2$ .

Exercise 3.8.4: Find a fundamental matrix solution for the system $x′1 = 7x1 + 4x2 + 12x3$ , $′ x2 = x1 + 2x 2 + x3$ , $′ x3 = -3x 1 - 2x2 - 5x3$ . Then find the solution that satisfies $[ ] 0 ⃗x = -12$ .

Exercise 3.8.5: Compute the matrix exponential $eA$ for $A = [1 2] 0 1$ .

Exercise 3.8.6: Suppose $AB = BA$ (matrices commute). Show that $eA+B = eAeB$ .

Exercise 3.8.7: Use exercise 3.8.6 to show that $(eA)-1 = e-A$ . In particular this means that $eA$ is invertible even if $A$ is not.

Exercise 3.8.8: Suppose $A$ is a matrix with eigenvalues $- 1$ , 1, and corresponding eigenvectors $[] 11$ , $[] 01$ . a) Find matrix $A$ with these properties. b) Find the fundamental matrix solution to $⃗x′ = A⃗x$ . c) Solve the system in with initial conditions $[ ] x⃗(0 ) = 23$ .

Exercise 3.8.9: Suppose that $A$ is an $n × n$ matrix with a repeated eigenvalue $λ$ of multiplicity $n$ . Suppose that there are $n$ linearly independent eigenvectors. Show that the matrix is diagonal, in particular $A = λI$ . Hint: Use diagonalization and the fact that the identity matrix commutes with every other matrix.

3.9 Nonhomogeneous systems

Note: 3 lectures (may have to skip a little), somewhat different from §5.6 in EP

3.9.1 First order constant coefficient

Integrating factor

Let us first focus on the nonhomogeneous first order equation

′ ⃗ ⃗x (t) = Ax⃗(t) + f (t),

where $A$ is a constant matrix. The first method we will look at is the integrating factor method. For simplicity we rewrite the equation as

⃗x′(t) + P⃗x(t) = f⃗(t),

where $P = -A$ . We multiply both sides of the equation by $etP$ (being mindful that we are dealing with matrices which may not commute) to obtain

etP⃗x′(t) + etPP ⃗x(t) = etPf⃗(t).

We notice that $tP tP Pe = e P$ . This fact follows by writing down the series definition of $tP e$ ,

( 1 ) 1
PetP = P I + I + tP +-(tP)2 + ⋅⋅⋅ = P + tP 2 +--t2P 3 + ⋅⋅⋅ =
2 2 ( )
1- 2 tP
= I + I + tP + 2 (tP) + ⋅⋅⋅ P = Pe .

We have already seen that $( ) d- etP = PetP dt$ . Hence,

-d( tP ) tP⃗ dt e ⃗x(t) = e f(t).

We can now integrate. That is, we integrate each component of the vector separately

∫ etP⃗x(t) = etP ⃗f(t) dt + ⃗c.

Recall from exercise 3.8.7 that $(etP)-1 = e-tP$ . Therefore, we obtain

∫ ⃗x(t) = e-tP etPf⃗(t) dt + e-tP⃗c.

Perhaps it is better understood as a definite integral. In this case it will be easy to also solve for the initial conditions as well. Suppose we have the equation with initial conditions

⃗x′(t) + P ⃗x(t) = ⃗f(t), ⃗x(0) = ⃗b.

The solution can then be written as

|----------∫---------------------| | -tP t sP ⃗ -tP⃗ | |⃗x(t) = e e f (s) ds + e b. | -------------0--------------------

(3.5)

Again, the integration means that each component of the vector $esPf⃗(s)$ is integrated separately. It is not hard to see that (3.5) really does satisfy the initial condition $⃗x(0) = ⃗b$ .

∫ 0 ⃗x(0) = e-0P esP ⃗f(s) ds + e-0P⃗b = I⃗b = ⃗b. 0

Example 3.9.1: Suppose that we have the system

x′+ 5x1 - 3x2 = et, 1′ x2 + 3x1 - x2 = 0,

with initial conditions $x1(0) = 1,x2(0 ) = 0$ .

Let us write the system as

[ ] [ t] [ ] ⃗x ′ + 5 -3 ⃗x = e , ⃗x(0) = 1 . 3 -1 0 0

We have previously computed $etP$ for $P = [5 -3] 3 -1$ . We immediately also have $e-tP$ , simply by negating $t$ .

[(1 + 3t)e2t - 3te2t ] [(1 - 3t)e-2t 3te-2t ] etP = 2t 2t , e-tP = -2t -2t . 3te (1 - 3t)e -3te (1 + 3t)e

Instead of computing the whole formula at once. Let us do it in stages. First

∫ t ∫ t[ 2s 2s ] [ s] esP⃗f(s) ds = (1 + 3s2)se - 3se 2s e ds 0 0 3se (1 - 3s)e 0 ∫ t[(1 + 3s)e3s] = 3s ds [0 3se] te3t = (3t-1)e3t+1-. 3

Then

∫ t ⃗x(t) = e-tP esP⃗f(s) ds + e-tP⃗b 0 [ -2t -2t ] [ te3t ] [ -2t -2t ] [ ] = (1 - 3t)-e2t 3te -2t (3t-1)e3t+1 + (1 - 3t)e-2t 3te -2t 1 - 3te (1 + 3t)e ----3--- - 3te (1 + 3t)e 0 [ te -2t ] [ -2t] = et (1 ) -2t + (1 - 3t)-e2t - 3 + 3 + t e -3te [ (1 - 2t)e-2t ] = et (1 ) -2t . - 3 + 3 - 2t e

Phew!

Let us check that this really works.

x′ + 5x1 - 3x2 = (4te-2t - 4e -2t) + 5(1 - 2t)e-2t + et - (1 - 6t)e-2t = et. 1

Similarly (exercise) $′ x2 + 3x 1 - x2 = 0$ . The initial conditions are also satisfied as well (exercise).

For systems, the integrating factor method only works if $P$ does not depend on $t$ , that is, $P$ is constant. The problem is that in general

d ∫ ∫ --e P(t)dt ⇔ P(t)e P(t)dt, dt

because matrices generally do not commute.

Eigenvector decomposition

For the next method, we note that the eigenvectors of a matrix give the directions in which the matrix acts like a scalar. If we solve our system along these directions these solutions would be simpler as we can treat the matrix as a scalar. We can put those solutions together to get the general solution.

Take the equation

⃗x′(t) = A ⃗x(t) +f⃗(t).

(3.6)

Assume that $A$ has $n$ linearly independent eigenvectors $⃗v ,...,⃗v 1 n$ . Let us write

⃗x (t) = ⃗v1ξ1(t) + ⃗v2ξ2(t) + ⋅⋅⋅ + ⃗vnξn(t).

(3.7)

That is, we wish to write our solution as a linear combination of the eigenvectors of $A$ . If we can solve for the scalar functions $ξ1$ through $ξn$ we have our solution $⃗x$ . Let us decompose $⃗ f$ in terms of the eigenvectors as well. Write

f⃗(t) = ⃗v g (t) + ⃗v g (t) + ⋅⋅⋅ + ⃗v g (t). 1 1 2 2 n n

(3.8)

That is, we wish to find $g1$ through $gn$ that satisfy (3.8). We note that since all the eigenvectors of $A$ are independent, the matrix $E = [⃗v1 ⃗v2 ⋅⋅⋅ ⃗vn]$ is invertible. We see that (3.8) can be written as $⃗f = E⃗g$ , where the components of $⃗g$ are the functions $g1$ through $gn$ . Then $⃗g = E -1⃗f$ . Hence it is always possible to find $⃗g$ when there are $n$ linearly independent eigenvectors.

We plug (3.7) into (3.6), and note that $A⃗vk = λk⃗vk$ .

′ ′ ′ ′ ⃗x = ⃗v1ξ(1 + ⃗v2ξ2 + ⋅⋅⋅ +⃗vnξn ) = A ⃗v1ξ1 + ⃗v2ξ2 + ⋅⋅⋅ + ⃗vnξn + ⃗v1g1 +⃗v2g2 + ⋅⋅⋅ + ⃗vngn = A⃗v1ξ1 + A⃗v2ξ2 + ⋅⋅⋅ + A ⃗vnξn + ⃗v1g1 + ⃗v2g2 + ⋅⋅⋅ + ⃗vngn = ⃗v1λ1ξ1 + ⃗v2λ 2ξ2 + ⋅⋅⋅ + ⃗vnλnξn + ⃗v1g1 + ⃗v2g2 + ⋅⋅⋅ +⃗vngn = ⃗v (λ ξ + g ) + ⃗v (λ ξ + g ) + ⋅⋅⋅ +⃗v (λ ξ + g ). 1 1 1 1 2 2 2 2 n n n n

If we identify the coefficients of the vectors $⃗v 1$ through $⃗v n$ we get the equations

ξ′= λ 1ξ1 + g1, 1′ ξ2 = λ 2ξ2 + g2, .. . ξ′n = λnξn + gn.

Each one of these equations is independent of the others. They are all linear first order equations and can easily be solved by the standard integrating factor method for single equations. That is, for example for the $th k$ equation we write

ξ′(t) - λkξk(t) = gk(t). k

We use the integrating factor $-λkt e$ to find that

[ ] d--ξ (t)e-λkt = e-λktg (t). dx k k

Now we integrate and solve for $ξk$ to get

∫ ξk(t) = eλkt e-λktgk(t) dt + Ckeλkt.

Note that if you are looking for just any particular solution, you could set $Ck$ to be zero. If we leave these constants in, we will get the general solution. Write $⃗x(t) = ⃗v1ξ1(t) + ⃗v2ξ2(t) + ⋅⋅⋅ +⃗vnξn(t)$ , and we are done.

Again, as always, it is perhaps better to write these integrals as definite integrals. Suppose that we have an initial condition $⃗x(0) = ⃗b$ . We take $⃗c = E-1⃗b$ and note $⃗b = ⃗v1a1 + ⋅⋅⋅ + ⃗vnan$ , just like before. Then if we write

|-----------∫--t---------------------| | ξ (t) = eλkt e-λksg (s) dt + a eλkt, | | k 0 k k | -------------------------------------

we will actually get the particular solution $⃗x(t) = ⃗v1ξ1(t) + ⃗v2ξ2(t) + ⋅⋅⋅ + ⃗vnξn(t)$ satisfying $⃗ ⃗x(0) = b$ , because $ξk(0) = ak$ .

Example 3.9.2: Let $[ ] A = 13 31$ . Solve $′ ⃗x = A⃗x + ⃗f$ where $[ t] f⃗(t) = 2e2t$ for $[ 3∕16] ⃗x(0) = -5∕16$ .

The eigenvalues of $A$ are $- 2$ and 4 and the corresponding eigenvectors are $[1] -1$ and $[1] 1$ respectively. This calculation is left as an exercise. We write down the matrix $E$ of the eigenvectors and compute its inverse (using the inverse formula for $2 × 2$ matrices)

[ ] [ ] 1 1 -1 1-1 - 1 E = -1 1 , E = 2 1 1 .

We are looking for a solution of the form $[ ] [ ] ⃗x = -11 ξ1 + 11 ξ2$ . We also wish to write $f⃗$ in terms of the eigenvectors. That is we wish to write $[ ] f⃗= 2et = [1]g + [1]g 2t -1 1 1 2$ . Thus

[ ] [ ] [ ] [ ] [ ] g1 -1 2et 1- 1 -1 2et et - t g2 = E 2t = 2 1 1 2t = et + t .

So $t g1 = e - t$ and $t g2 = e + t$ .

We further want to write $⃗x(0)$ in terms of the eigenvectors. That is, we wish to write $[ ] ⃗x(0) = 3∕16 = [1] a + [1]a -5∕16 -1 1 1 2$ . Hence

[ ] [ ] [ ] a1 -1 316 14 = E -5 = -1 . a2 16 16

So $a1 = 1∕4$ and $a2 = -1∕16$ . We plug our $⃗x$ into the equation and get that

[ ] [ ] [ ] [ ] [ ] [ ] 1 ξ′1 + 1 ξ2′= A 1 ξ1 + A 1 ξ2 + 1 g1 + 1 g2 -1 1 -1 1 - 1 1 [ 1 ] [1] [ 1 ] [1 ] = (-2ξ1) + 4ξ2 + (et - t) + (et - t). - 1 1 -1 1

We get the two equations

′ t 1- ξ1 = - 2ξ1 + e - t, where ξ1(0) = a1 = 4, -1 ξ′2 = 4 ξ2 + et + t, where ξ2(0) = a2 = ---. 16

We solve with integrating factor. Computation of the integral is left as an exercise to the student. Note that you will need integration by parts.

∫ et t 1 ξ1 = e-2t e2t(et - t) dt + C 1e-2t =----+ --+ C 1e-2t. 3 2 4

$C 1$ is the constant of integration. As $ξ1(0) = 1∕4$ then $1∕4 = 1∕3 + 1∕4 + C 1$ and hence $C 1 = -1∕3$ . Similarly

∫ 4t -4t t 4t et -t 1-- 4t ξ2 = e e (e + t) dt + C2e = - 3 - 4 - 16 + C2e .

As $ξ (0) = 1∕16 2$ we have that $-1∕16 = -1∕3 - 1∕16 + C 2$ and hence $C = 1∕3 2$ . The solution is

[ ] (et - e-2t 1 - 2t) [ ] (e4t - et 4t + 1 ) [ e4t-e-2t- 3-12t ] ⃗x(t) = 1 --------+ ------ + 1 -------- ------ = e-2t+3e4t+2e+t 164t-5 . -1 3 4 1 3 16 ---3----+ 16-

That is, $x1 = e4t-e-2t+ 3-12t 3 16$ and $x2 = e-2t+e4t+2et+ 4t-5 3 16$ .

Exercise 3.9.1: Check that $x 1$ and $x 2$ solve the problem. Check both that they satisfy the differential equation and that they satisfy the initial conditions.

Undetermined coefficients

The method of undetermined coefficients also still works. The only difference here is that we will have to take unknown vectors rather than just numbers. Same caveats apply to undetermined coefficients for systems as they do for single equations. This method does not always work, furthermore if the right hand side is complicated, you will have lots of variables to solve for. In this case you can think of each element of an unknown vector as an unknown number. So in system of 3 equations if you have say 4 unknown vectors (this would not be uncommon), then you already have 12 unknowns that you need to solve for. The method can turn into a lot of tedious work. As this method is essentially the same as it is for single equations, let us just do an example.

Example 3.9.3: Let $[ ] A = --12 01$ . Find a particular solution of $′ ⃗x = A ⃗x + f⃗$ where $[t] f⃗(t) = et$ .

Note that we can solve this system in an easier way (can you see how), but for the purposes of the example, let us use the eigenvalue method plus undetermined coefficients.

The eigenvalues of $A$ are $- 1$ and 1 and the corresponding eigenvectors are $[1] 1$ and $[0] 1$ respectively. Hence our complementary solution is

[1] [0] ⃗xc = α 1 e-t + α2 et, 1 1

for some arbitrary constants $α 1$ and $α 2$ .

We would want to guess a particular solution of

⃗x = ⃗aet + ⃗bt + ⃗c.

However, something of the form $⃗aet$ appears in the complementary solution. Because we do not yet know if the vector $⃗a$ is a multiple of $[] 01$ we do not know if a conflict arises. It is possible that no conflict arises, but to be safe we should also try $⃗btet$ . Here we find the crux of the difference for systems. You want to try both $t ⃗ae$ and $⃗ t bte$ in your solution, not just $⃗ t bte$ . Therefore, we try

⃗x = ⃗aet + ⃗btet + ⃗ct + ⃗d.

Thus we have 8 unknowns. We write $⃗a = [a1] a2$ , $⃗b = [b1] b2$ , $⃗c = [c1] c2$ , and $⃗d = [d1] d2$ , We have to plug this into the equation. First let us compute $′ ⃗x$ .

⃗x′ = (⃗a + ⃗b)et +⃗btet + ⃗c.

Now $′ ⃗x$ must equal $A ⃗x + f⃗$ so

⃗ t ⃗ t ⃗ ⃗
A⃗x + f = A⃗ae + A bte +[A⃗ct + A d]+f = [ ] [ ] [ ] [ ]
- a1 t -b1 t -c1 - d1 et
= - 2a + a e + -2b + b te + - 2c + c t + - 2d + d + t .
1 2 1 2 1 2 1 2

Now we identify the coefficients of $t e$ , $t te$ , $t$ and any constants.

a 1 + b1 = -a1 + 1, a 2 + b2 = -2a1 + a2, b = -b , 1 1 b2 = -2b1 + b2, 0 = -c , 1 0 = -2c1 + c2 + 1, c1 = -d1, c2 = -2d1 + d2.

We could write this is an $8 × 9$ augmented matrix and start row reduction, but it is easier to just do this in an ad hoc manner. Immediately we see that $b1 = 0$ , $c1 = 0$ , $d1 = 0$ . Plugging these back in we get that $c = -1 2$ and $d = -1 2$ . The remaining equations that tell us something are

a = - a + 1, 1 1 a2 + b2 = - 2a1 + a 2.

So $a1 = 1 2$ and $b2 = - 1$ . $a2$ can be arbitrary and still satisfy the equation. We are looking for just a single solution so presumably the simplest one is when $a = 0 2$ . Therefore,

[1] [ ] [ ] [ ] [ 1 t ] ⃗x = ⃗aet + ⃗btet + ⃗ct + ⃗d = 2 et + 0 tet + 0 t + 0 = t2e . 0 -1 - 1 -1 -te - t - 1

That is, $x1 = 1et 2$ , $x 2 = - tet - t - 1$ . You would add this to the complementary solution to get the general solution of the problem. Notice also that both $t ⃗ae$ and $⃗ t bte$ really was needed.

Exercise 3.9.2: Check that $x 1$ and $x 2$ solve the problem. Also try setting $a = 1 2$ and again check these solutions. What is the difference between the two solutions we can obtain in this way?

As you can see, other than the handling of conflicts, undetermined coefficients works exactly the same as it did for single equations. However, the computations can get out of hand pretty quickly for systems. The equation we had done was very simple.

3.9.2 First order variable coefficient

Just as for a single equation, there is the method of variation of parameters. In fact for constant coefficient systems, this is essentially the same thing as the integrating factor method we discussed earlier. However this method will work for any linear system, even if it is not constant coefficient, provided you have somehow solved the associated homogeneous problem.

Suppose we have the equation

′ ⃗x = A (t)⃗x + f⃗(t).

(3.9)

Further, suppose that you have solved the associated homogeneous equation $′ ⃗x = A(t)⃗x$ and found the fundamental matrix solution $X(t)$ . The general solution to the associated homogeneous equation is $X (t)⃗c$ for a constant vector $⃗c$ . Just like for variation of parameters for single equation we try the solution to the nonhomogeneous equation of the form

⃗xp = X (t)⃗u(t),

where $⃗u (t)$ is a vector valued function instead of a constant. Now substitute into (3.9) to obtain

⃗xp′(t) = X ′(t)⃗u(t) + X(t)⃗u ′(t) = A (t)X(t)⃗u (t) + f⃗(t).

But $X$ is the fundamental matrix solution to the homogeneous problem so $′ X (t) = A (t)X(t)$ , and thus

X ′(t)⃗u (t) + X (t)⃗u′(t) = X′(t)⃗u(t) + ⃗f(t).

Hence $X(t)⃗u ′(t) = ⃗f(t)$ . If we compute $[X(t)]-1$ , then $⃗u′(t) = [X (t)]-1 ⃗f(t)$ . Now integrate to obtain $⃗u$ and we have the particular solution $⃗xp = X(t)⃗u (t)$ . Let us write this as a formula

|----------∫-----------------| | ⃗xp = X (t) [X (t)]-1 ⃗f(t) dt. ------------------------------

Note that if $A$ is constant and you let $tA X(t) = e$ , then $-1 -tA [X(t)] = e$ and hence we get a solution $∫ ⃗xp = etA e-tAf⃗(t) dt$ which is precisely what we got using the integrating factor method.

Example 3.9.4: Find a particular solution to

′ 1 [t -1] [t] ⃗x = -2---- ⃗x + (t2 + 1 ). t + 1 1 t 1

(3.10)

Here $[ ] A = t21+1 t1 -t1$ is most definitely not constant. Perhaps by a lucky guess, you find that $[ ] X = 1t -1t$ solves $X ′(t) = A(t)X (t)$ . Once we know the complementary solution we can easily find a solution to (3.10). First we find

[ ] [X(t)]-1 = --1--- 1 t . t2 + 1 - t 1

Next we know a particular solution to (3.10) is

∫ ⃗xp = X(t) [X(t)]-1f⃗(t) dt [ ]∫ [ ] [ ] 1 - t --1--- 1 t t 2 = t 1 t2 + 1 - t 1 1 (t + 1) dt [ ]∫ [ ] 1 - t 2t = t 1 - t2 + 1 dt [ ] [ 2 ] = 1 - t 1 t3 t 1 -3t + t [ 1t4 ] = 2 33 . 3t + t

Adding the complementary solution we have that the general solution to (3.10).

[ ] [ ] [ 1 4 ] [ 1 4 ] ⃗x = 1 - t c1 + 2 33t = c1 - c2t + 3t 23 . t 1 c2 3t + t c2 + (c1 + 1)t + 3t

Exercise 3.9.3: Check that $x = 1t4 1 3$ and $x = 2t3 + t 2 3$ really solve (3.10).

In the variation of parameters, just like in the integrating factor method we can obtain the general solution by adding in constants of integration. That is, we will add $X (t)⃗c$ for a vector of arbitrary constants. But that is precisely the complementary solution.

3.9.3 Second order constant coefficients

Undetermined coefficients

We have already previously did a simple example of the method of undetermined coefficients for second order systems in § 3.6. This method is essentially the same as undetermined coefficients for first order systems. There are some simplifications that you can make, as we did in § 3.6. Let the equation be

′′ ⃗x = A ⃗x + F⃗(t),

where $A$ is a constant matrix. If $⃗F(t)$ is of the form $⃗F 0cosωt$ , then you can try a solution of the form

⃗xp = ⃗cco sωt,

and you do not need to introduce sines.

If the $⃗F$ is a sum of cosines, you note that we still have the superposition principle, so if $F⃗(t) = ⃗F0 cosω 0t + ⃗F1cos ω1t$ , you could try $⃗acos ω0t$ for the problem $⃗x′′ = A ⃗x + F⃗0cos ω0t$ , and you would try $⃗ bcos ω1t$ for the problem $′′ ⃗ ⃗x = A ⃗x + F 0cosω 1t$ . Then sum the solutions.

However, if there is duplication with the complementary solution, or the equation is of the form $′′ ′ ⃗x = A⃗x + B ⃗x + F⃗(t)$ , then you need to do the same thing as you do for first order systems.

Actually you will never go wrong with putting in more terms than needed into your guess. You will just find that the extra coefficients will turn out to be zero. But it is useful to save some time and effort.

Eigenvector decomposition

If we have the system

⃗x′′ = A ⃗x + F⃗(t),

we can do eigenvector decomposition, just like for first order systems.

Let $λ1$ , …, $λn$ be the eigenvalues and $⃗v1$ , …, $⃗vn$ be the eigenvectors. Again form the matrix $E = [⃗v1 ⋅⋅⋅⃗vn]$ . Write

⃗x (t) = ⃗v1ξ1(t) + ⃗v2ξ2(t) + ⋅⋅⋅ + ⃗vnξn(t).

Decompose $⃗F$ in terms of the eigenvectors

⃗F (t) = ⃗v1g1(t) + ⃗v2g2(t) + ⋅⋅⋅ +⃗vngn(t).

And again $⃗g = E -1⃗F$ .

Now plug in and doing the same thing as before

′′ ′′ ′′ ′′ ⃗x = ⃗v1ξ(1 + ⃗v2ξ2 + ⋅⋅⋅ +⃗vnξn ) = A ⃗v1ξ1 + ⃗v2ξ2 + ⋅⋅⋅ + ⃗vnξn + ⃗v1g1 + ⃗v2g2 + ⋅⋅⋅ + ⃗vngn = A ⃗v ξ + A⃗vξ + ⋅⋅⋅ + A⃗v ξ + ⃗v g + ⃗v g + ⋅⋅⋅ + ⃗v g 1 1 22 n n 1 1 2 2 n n = ⃗v1λ 1ξ1 + ⃗v2λ2ξ2 + ⋅⋅⋅ + ⃗vnλnξn + ⃗v1g1 +⃗v2g 2 + ⋅⋅⋅ + ⃗vngn = ⃗v1(λ1ξ1 + g 1) + ⃗v2(λ 2ξ2 + g2) + ⋅⋅⋅ + ⃗vn(λnξn + gn).

Identify the coefficients of the eigenvectors to get the equations

′′ ξ1 = λ1ξ1 + g1, ξ′′= λ2ξ2 + g2, 2 ... ′′ ξn = λnξn + gn.

Each one of these equations is independent of the others. Now solve each one of these using the methods of chapter 2. Now write $⃗x(t) = ⃗v1ξ1(t) + ⋅⋅⋅ + ⃗vnξn(t)$ , and we are done; we have a particular solution. If you have found the general solution for $ξ1$ through $ξn$ , then again $⃗x(t) = ⃗v1ξ1(t) + ⋅⋅⋅ + ⃗vnξn(t)$ is the general solution.

Example 3.9.5: Let us do the example from § 3.6 using this method. The equation is

[ ] [ ] ′′ -3 1 0 ⃗x = 2 -2 ⃗x + 2 cos3t.

The eigenvalues were $- 1$ and $- 4$ , with eigenvectors $[] 12$ and $[ ] 1-1$ . Therefore $[ ] E = 12 1-1$ and $[ ] E -1 = 13 12 -11$ . Therefore,

[ ] [ ] [ ] [ 2 ] g 1 -1⃗ 1- 1 1 0 3 co s3t g 2 = E F (t) = 3 2 -1 2cos 3t = -2-cos3t . 3

So after the whole song and dance of plugging in, the equations we get are

2 ξ′′1 = -ξ1 + -co s3t, 3 ξ′′ = -4ξ - 2-cos3t. 2 2 3

For each we can try the method of undetermined coefficients and try $C 1cos 3t$ for the first equation and $C 2cos3t$ for the second equation. We plug in to get

2 - 9C 1cos3t = -C 1cos3t + --cos3t, 3 - 9C 2cos3t = -4C 2cos 3t - 2-cos3t. 3

Each of these equations we solve separately: we get $2 - 9C 1 = -C 1 + ∕3$ and $2 - 9C 2 = - 4C2 - ∕3$ . And hence $C 1 = -1∕12$ and $C2 = 2∕15$ . So our particular solution is

[ ] ( ) [ ] ( ) [1 ] ⃗x = 1 -1-cos3t + 1 -2-cos3t = ∕20 co s3t. 2 12 -1 15 -3∕10

This solution matches what we got previously in § 3.6.

3.9.4 Exercises

Exercise 3.9.4: Find a particular solution to $′ x = x + 2y + 2t$ , $′ y = 3x + 2y - 4$ , a) using integrating factor method, b) using eigenvector decomposition, c) using undetermined coefficients.

Exercise 3.9.5: Find the general solution to $x′ = 4x + y - 1$ , $y′ = x + 4y - et$ , a) using integrating factor method, b) using eigenvector decomposition, c) using undetermined coefficients.

Exercise 3.9.6: Find the general solution to $′′ x1 = -6x 1 + 3x2 + cost$ , $′′ x2 = 2x1 - 7x2 + 3 cost$ , a) using eigenvector decomposition, b) using undetermined coefficients.

Exercise 3.9.7: Find the general solution to $x′′1 = -6x 1 + 3x2 + co s2t$ , $x′2′= 2x1 - 7x2 + 3 cos2t$ , a) using eigenvector decomposition, b) using undetermined coefficients.

Exercise 3.9.8: Take the equation

[1 ] [ 2] ⃗x′ = t -11 ⃗x + t . 1 t - t

a) Check that

[ ] [ ] tsin t tcost ⃗xc = c1 -tcos t + c2 tsin t

is the complementary solution. b) Use variation of parameters to find a particular solution.

Chapter 3Systems of ODEs

3.1 Introduction to systems of ODEs

3.1.1 Exercises

3.2 Matrices and linear systems

3.2.1 Matrices and vectors

3.2.2 Matrix multiplication

3.2.3 The determinant

3.2.4 Solving linear systems

3.2.5 Computing the inverse

3.2.6 Exercises

3.3 Linear systems of ODEs

3.3.1 Exercises

3.4 Eigenvalue method

3.4.1 Eigenvalues and eigenvectors of a matrix

3.4.2 The eigenvalue method with distinct real eigenvalues

3.4.3 Complex eigenvalues

3.4.4 Exercises

3.5 Two dimensional systems and their vector fields

3.5.1 Exercises

3.6 Second order systems and applications

3.6.1 Undamped mass spring systems

3.6.2 Examples

3.6.3 Forced oscillations

3.6.4 Exercises

3.7 Multiple eigenvalues

3.7.1 Geometric multiplicity

3.7.2 Defective eigenvalues

3.7.3 Exercises

3.8 Matrix exponentials

3.8.1 Definition

3.8.2 Simple cases

3.8.3 General matrices

3.8.4 Fundamental matrix solutions

3.8.5 Approximations

3.8.6 Exercises

3.9 Nonhomogeneous systems

3.9.1 First order constant coefficient

Integrating factor

Eigenvector decomposition

Undetermined coefficients

3.9.2 First order variable coefficient

3.9.3 Second order constant coefficients

Undetermined coefficients

Eigenvector decomposition

3.9.4 Exercises

Chapter 3
Systems of ODEs