4. Systems of ODEs

§4.1 Introduction

    In the Assignment 2 we used PPLANE to explore 2x2 systems of ODEs.  In this lab we further investigate system of equations and focus on linear ones, the advantages of specialization being that the techniques of eigenvalues and eigenvectors you will study here are extremely effective even for linear systems of large size.   A homogeneous 2x2 system of linear ODE's has the form

(1) x' = ax+ by

y' = cx+ dy

where a, b, c, and d are constants.

    To solve systems of linear differential equations, it is often helpful to "rephrase" the problem in matrix notation.  The above system can be expressed as

v' = Av

where v is the column vector   and A is the matrix .  Such systems can be solved using the eigenvalues and eigenvectors of the matrix A as we shall see below.  To begin with, let's look briefly at how matrices are entered and manipulated in MATLAB.

§4.2 Matrices in MATLAB

    To enter a matrix into MATLAB, we use square brackets to begin and end the contents of the matrix and semicolons to separate the rows.  The components of a single row are separated by commas.  For example, the command

>> A = [1, 2; 3, 4]

will result in the assignment of a matrix to the variable A:
 

A =

    1    2

    3    4

We can enter a column vector by thinking of it as an m×1 matrix, so the command

>> b = [1; -1]

will result in the 2×1 column vector:
 

b =

    1
   -1

 

    There are many properties of matrices that MATLAB will calculate through simple commands.  Most of this material is covered in Math 20F, but a few basic properties can be helpful to us in solving systems of differential equations. 

An eigenvalue s is related to its eigenvector b by the equation

            Ab = sb.

MATLAB can be used to find the eigenvalues and eigenvectors of a matrix using the eig command.  The command by itself, eig(A), will return a column vector with the eigenvalues of A as its components; for example, for our matrix A above, we get the following output:

>> eig(A)

ans =

        -0.3723
        5.3723

If we also want MATLAB to calculate the eigenvectors of A, we need to specify two output variables.  The command and its output are below:

>> [eigvec, eigval] = eig(A)

eigvec =

        -0.8246     -0.4160
         0.5658     -0.9094


eigval =

        -0.3723     0
         0          5.3723

With this command, we get two matrices as output.  The first matrix, called eigvec, has the eigenvectors of A as its columns.  The second matrix, eigval, has the eigenvalues of A on its main diagonal and zeros everywhere else.

Let B be the matrix

B =

   

Exercise 4.1

Include input and output.

 (a)  Define the matrix B in MATLAB with the values above.  Copy and paste the input and output from your command into your Word document.

Include input and output.

 (b)  Use the MATLAB commands to find the eigenvalues and eigenvectors for the matrix B.  Copy and paste the input and output from your command into your Word document.

§4.3 Solving a Homogeneous System of Linear ODEs

    We are now ready to solve systems of equations such as (1) above.

    As you have probably already seen in class, when the matrix A = has two distinct eigenvalues, the general solution to the system (1) is

(2)

where s1 and s2 are the eigenvalues of A, and b1 and b2 are corresponding eigenvectors, and c1 and c2 are constants, which are determined by the initial conditions.  It is important to note that the eigenvalues can be real or complex.  In the remainder of this lab, we will examine how the eigenvalues of a linear system determine the behavior of the solution curves.

Real Eigenvalues

    Let us begin by examining a what happens when the eigenvalues are real.

Example 4.1

    Now that we have the general formula for the solutions to a 2x2 system of differential equations, let's see it used in an example. Consider the system of differential equations

   
(3) x' = -2.1x + 1.2y

y' = 0.8x -3.4y

When the system is put into the form v' = Av, the matrix is given by A = .

Enter this matrix into MATLAB with the command:

>> A = [-2.1, 1.2; 0.8, -3.4];

We can get the eigenvalues and eigenvectors of A by typing

>> [eigvec, eigval] = eig(A)

eigvec =

    0.9159   -0.5492
    0.4013    0.8357

eigval =

   -1.5742         0
         0   -3.9258

Using the eigenvalues and eigenvectors listed above, we find the general solution



  Note that we can tell quite a bit about the solution just by looking at it qualitatively. Since both eigenvalues are negative, and since they constitute the rate constants of the exponentials, the solution will tend to 0 as t gets larger regardless of which values we choose for c1 and c2, .

We can also get a visual representation of the solutions by using PPLANE (revisit Assignment 2 for more information about PPLANE). On the PPLANE display window, choose Options -> Solution direction -> Forward. This will show us how the solutions behave as t gets large. Clicking a few places on the window, we get something that looks like:

It is worth noting that the solutions we got from PPLANE all tend to 0, just as we predicted by looking at the eigenvalues.

    Now consider the system of differential equations

   
(4) x' = x + 3y

y' = -x - 8y

 

Exercise 4.2
Do work by hand.
 (a)  When the system (4) above is put into the form v' = Av, what is the matrix A?  Leave space in your Word document to write this in.
Include input and output.

(b)  Enter the matrix A into MATLAB.  Use MATLAB to find the characteristic roots (eigenvalues) and characteristic vectors (eigenvectors) of your matrix A.  Copy and paste all input and output from your command into your Word document.
Do work by hand. (c)  Use (2) above, and the results from part (b), to write the general solution of our system (4).  Leave space in your Word document for this.  What happens to the system as t gets large?
Do work by hand. (d)  Use PPLANE to create a plot of the solutions. As in the example, choose Options -> Solution direction -> Forward and click on the plot several times. Include this plot in your writeup.
Does this match your answer to (c)?

§4.3 Complex Eigenvalues

   In our 2x2 systems thus far, the eigenvalues and eigenvectors have always been real numbers. However, it is entirely possible that both the eigenvalues and the eigenvectors of a 2x2 matrix will have complex entries. As long as the eigenvalues are distinct, we will still have a general solution of the form given above in (2).  In order to build up an understanding of this equation, we will take a brief deviation for a quick review of complex exponentials.

Review of Complex Exponentials and Complex Conjugates

    If our system has complex eigenvalues, then formula (2) tells us we will have to deal with complex exponentials, which we will now briefly review.  Recall that if we can write.
Splitting up the exponential like this allows us  to see how the real and imaginary parts of s effect the exponential.  If r < 0, the function will damp out as t gets larger.   If, on the other hand, r >0, will cause the exponential to explode as t gets larger.  The purely complex exponential, , has absolute value one; as t changes, it will move in a circle when plotted on the complex plane.  The bigger , the faster the the function circles.  Thus the complex exponential is an example of damped or exploding oscillatory motion.

   We shall be interested only in real valued functions v(t), so let's review some complex vs real facts. Recall that if is a complex number, then its complex conjugate is given by .  The real part r of s is half of s plus s conjugate. Indeed, v(t) in (2) is real if and only if the second term in (2) is the complex conjugate of the first. Thus, if we have 2x2 system with complex eigenvalues and , with corresponding eigenvectors and , then we have that (2) is given by

See page 404 of Boyce & DiPrima for further details. Note that the "Re" here takes the real part of the complex number, and "Im" takes the imaginary part.  If we look at only the first term of the above, we see that this can be written as

where we have written out the eigenvector as a column vector. Let's get some practice seeing what these individual solutions look like.

Example 4.2

Consider the complex exponential given by .  The graph of this looks like:


It is important to realize the vector itself can be complex.  Consider for instance the function  which has graph


    Complex Eigenvalues and Differential Equations

    Typically in practical problems things such as data and measurements consist of real numbers rather than complex numbers. However, when we solve for the eigenvalues and eigenvectors of A, we may get complex values, which may seem counterintuitive. The following theorem makes sense of this:

Theorem 4.1 (see Boyce and DiPrima section 7.6):

If the entries of A are real numbers and if the initial condition is a vector with real entries, then the solution v(t) is a vector with real entries.

   Let us see what algebraic structure this theorem actually says. Suppose that the matrix A has an eigenvector b with corresponding eigenvalue s, so that Ab = sb. One consequence of the theorem is that the complex conjugate of b is also an eigenvector, and its corresponding eigenvalue is the complex conjugate of s. This behavior ensures that even when the eigenvalues and eigenvectors have complex numbers in them, the solutions will still be real valued.

   The conclusion of the theorem can be very deceptive! Many students try to conclude that because the solutions will be real, we can ignore complex numbers. This is not true! As we will see later, complex eigenvalues lead to solutions which "spiral".  Indeed, complex numbers arise any time we have a system of differential equations describing ocsillatory motion.

Example 4.3

Now that we've built up our understanding of complex exponentials, let's do an example applying them to a 2x2 system. Consider the system of differential equations

   
(5) x' = -0.5x + 0.9y

y' = -1.2x -1.7y

Putting the system in the form v' = Av, we get the matrix A = .

We can use MATLAB to calculate the eigenvalues and eigenvectors as before. The output is

>> [eigvec, eigval] = eig(A)

eigvec =

  -0.3780 - 0.5345i  -0.3780 + 0.5345i
   0.7559             0.7559          

eigval =

  -1.1000 + 0.8485i        0          
        0            -1.1000 - 0.8485i

We have stated that the general solution of (2) is still valid for complex eigenvalues and eigenvectors, so we can find the general solution

 .

 Next plot the solutions using PPLANE and choose Options -> Solution direction -> Forward. Clicking several places on the plot, we get a plot that looks like:

See if you can identify which feature of the eigenvalues and/or eigenvectors dictates that the solutions tend to zero as t gets large.

   Not all solutions to 2x2 systems with complex eigenvalues and eigenvectors look like the plot in the last example. For further practice, take the system

   
(6) x' = 2.7x  - y

y' = 4.2x + 3.5y


Exercise 4.3
Include input and output.

(a)  Put the system above in the form v' = Av and enter the matrix A into MATLAB.  Use MATLAB to find eigenvalues and eigenvectors of A.  Copy and paste all input and output from your command into your Word document.
Do work by hand. (c)  Write the general solution of the system (6) by hand.  Leave space in your Word document for this.
Do work by hand.Include input and output. (c)  Use PPLANE to create a plot of the solutions. As in the example, choose Options -> Solution direction -> Forward and click on the plot several times. Include this plot in your writeup. What feature of the eigenvalues or eigenvectors causes the solutions to tend to infinity as t gets large?


Systems Larger Than 2x2

    So far, our discussion has focused only on 2x2 systems.  While we are unable to plot larger systems, our analysis of the eigenvalues still allows us to analyze this behavior.   For instance, a 4x4 matrix A with real entries typically has 4 eigenvalues s1, s2, s3, and s4, with corresponding eigenvectors b1, b2, b3, b4.  These correspond to solutions of the form:



The solution is real valued and its form depends on how many of the eigenvalues are real.

For example, if  s1 and  s2 are real but s3 and s4 are complex, then s3 and s4 are complex conjugates and all real solutions have the form
(7)

Another possibility is that none of the eigenvalues are real, so that we have two pairs of complex conjugates.  Suppose s1 is conjugate to s2 and s3  is conjugate to s4 .  Then our solution will have the form
(8)

What about s1 real and s2, s3, s4 not real? Is this possible? We leave this to your thoughts.


§4.4 Stability
One of the most important properties one worries about when dealing with situations govered by differential equations is whether or not the system is stable.  Roughly speaking, if you were to buy a piece of electronics which is not stable, then it will blow up when you plug it into the wall or produce large and unwanted oscillations.   Stability is not limited to electronics, however.  Economic models also can become unstable, as was seen in the financial meltdown of September 2008.  Of the many possibilities, we picked as an exercise in this class a mechanical system (an airplanes sensitivity to turbulence) because the phenomenon is intuitive to us all. Even designers of algorithms worry about the "stability of the algorithm"; convergence of an algorithm from a mathematical point of view amounts to a type of stability.

In terms of linear systems of differential equations, a linear system Ax = x' is called stable if any solution x(t) stays bounded as .  A linear system is called asymptotically stable if any solution as .  Asymptotical stability is stronger than stability; that is, stability is implied by asymptotic stablity.  In this section, we will examine stability for linear systems and how see how we can use the eigenvalues of a system to determine if it is stable, asymptotically stable, or unstable.  We will start with 2x2 systems, and then indicate how our analysis scales to systems of any size.

Example 4.4

    In Example 4.1, we saw that the system given by (3) had eigenvalues -1.5742 and -3.9258.  The analysis given there and the PPLANE plot of the direction field both showed that all the solutions tended to 0 as t got large.  Thus (3) is an example of a system which is asymptotically stable.  In fact, the eigenvalues of a system are both negative if and only if the system is asymptotically stable.  Contrast this with Exercise 4.2, in which you saw how the solutions to the system given in (4) got larger as .  This is an example of an unstable system. Notice the eigenvalues for this system are 0.6533 and -7.6533.  Any system with even one positive eigenvalue will be unstable.

For an additional example, consider the system given by

   
(9) x' = 2x - 5y

y' = x - 2y

A PPLANE plot looks like



Notice that the plot does not approach zero, but remains bounded.  Thus this is an example of a stable system which is not asymptotically stable.   MATLAB tells us that the eigenvalues of this system are i and -i.  The key here is that there is no real component to the eigenvalues.  Can you see, in terms of (2), how this determines the behavior of the solutions?

The thing to take away from this is that the eigenvalues determine the stability or non-stability of a system.  There are two important things to note. First, we have not considered all possible cases for the eigenvalues.  Second, while PPLANE is helpful for 2x2 systems, we cannot plot larger systems.  Thus we must look at the solution form (2) and do our analysis directly.  The following exercise explores this idea.

Consider the 3x3 system

(10) x'= 1.25x - 0.97y + 4.6z

y' = -2.6x - 5.2y - 0.31z

z' = 1.18x - 10.3y + 1.12z

Exercise 4.4

Include input and output.

 (a)  Define the matrix A in MATLAB with the values above, using A = [1.25, -.97, 4.6; -2.6, -5.2, -.31; 1.18, -10.3, 1.12], and calculate the eigenvalues of A. Copy and paste the input and output from MATLAB.

Include input and output.

 (b)  Using your answer from part (a), is the system in (10) stable? Justify your answer.

§4.5 Stability in Airplane Controls

    A fundamental property of all linear dynamical systems (meaning most objects which move) is that they have resonances. These are associated with eigenvaues and eigenvectors of the coeffiecient matrix of the system. While we could illustrate this with fluids in pipes, electrical circuits, signals in the air structures in earthquakes, etc., we have chosen to illustrate this with airplane since it is easy to visualize and something with which you are familiar.

We will consider a model used in the design of commercial aircraft. This is a differential equation which describes the effect of rate of change of rudder angle on the rate of change of yaw (the yaw is the side to side rotation of the aircraft; see the animation below).

The model in the next exercise is of a Boeing 747 at cruising speed assuming no control system is installed. It is a good thing we have a mathematical model, since such a plane could not fly for long, as you will soon see. The "yaw rate model" is a 4 x 4 system

(11) x'(t)   =   A x(t) + B u(t)

with a single function u as driving term. Here u is the rate of change of rudder angle and the 4 components of x are:

x1 = lateral (side to side) velocity
x2 = yaw rate
x3 = roll rate
x4 = pitch rate

         

The depiction above of yaw, roll and pitch respectively is courtesy of Wikipedia user ZeroOne, NASA and Wikipedia.

A is the 4 x 4 matrix

and B is the 4 x 1 matrix which is the transpose of (.01   -.175   .153   0). In our exercise we will not examine the effects of the driving term but instead just consider the homogeneous system x'= A x.

Exercise 4.5

Include input and output.

 (a)  What are the eigenvalues of the matrix A?

Include input and output.

 (b)  According to the mathematical definition is the system  x'=Ax   stable, asymptotically stable, or unstable?  Why?
Include input and output. (c)  One of the eigenvalues you obtained is very close to 0. Look at the corresponding eigenvector. Which component is biggest? Which phenomenon is this eigenvector most closely associated with: yaw, roll or pitch?

We think you might find the practical implications of the eigenvalue locations interesting, since 3 out of 4 of them have very small real parts (not strongly positive or strongly negative but at the transition):

(a) You found two real eigenvalues with one sufficiently negative that its effect damps out quickly; indeed, it has little effect on the performance of the airplane.

(b) The other real eigenvalue was close to zero and real. This corresponds to the plane being responsive to the pilots' commands, a property which you definitely want.  Note no oscillation is involved, since the eigenvalue is real.

(c) There is a pair of complex eigenvalues, say s3 and its complex conjugate, close to the imaginary axis. Even though their real part is negative, it is so small that their effect damps out slowly. They correspond to a basic tendency for the B747 to oscillate back-forth and roll at speed equal to Imaginary part of s3 radians per second. This is called the Dutch Roll Mode and can be very bad. A wind gust can "excite" the Dutch Roll Mode (which, if you have flown, you have almost certainly felt) and then another wind gust might well sustain (or increase) it.

(a)(b)(c) are often called resonant modes or eigenmodes of the airplane. For this airplane model there are 3 resonant modes.

The image above of the Dutch Roll Mode is courtesy of Wikipedia user Picascho and Wikipedia. 

Fortunately, this model is of the airframe alone, and not the airframe with the control system. Of course, engineers design (based on this mathematical model) and implement a control algorithm, called a yaw damper, which automatically moves the rudder back and forth and compensates for this phenomena. Large commercial airplanes require a yaw damper.

To see the Dutch Roll phenomenon, a 30 second video is Dutch Roll Video.  Keep in mind as you watch: the pilot is not touching the controls, but rather you are viewing a natural motion of the airplane.

Controlling the Dutch Roll

What is involved in making this plane reasonable to fly? At any time t, sensors tell us the state x(t) of the plane, and (roughly speaking) we can at that time design the rate of change of rudder angle u(t) to be anything we want.  Let us make it a simple function of the form

where F is a 1x4 matrix F = (F1, F2, F3, F4). Here is what the pilot tells the rudder she wants to occur and

Fx(t) =  F1 *x 1 + F2 *x 2 + F3 *x 3 + F4 *x 4

is what the plane's computer tells the rudder to do, with the aim of damping your planes bad resonant oscillations (which are often provoked by wind gusts). Our design amounts to choosing four real numbers F1, F2, F3, and F4. Putting this together with the airframe dynamical model given by (11) and we get

Stability of this system is completely determined by the eigenvalues of (A+BF). 

For the following exercise, recall matrices A and B are given right above exercise 4.5


Exercise 4.6

Include input and output.

 (a)  If F = (0 7 0 -1), what is the matrix BF?

Include input and output.

 (b)  If F = (0 5 0 -.1), what is the matrix BF?
Include input and output. (c)  If F = (0 5 0 -.1), what is the matrix A + BF? Are its eigenvalues the same or different from those of A?
Include input and output. (c)  Find F = (F1 F2 F3 F4) that increases the damping of the yaw oscillation and so that the plane still responds decently to pilot commands. More specifically: we want the complex eigenvalues to have real part less that -0.2 and that there is a real eigenvalue within 0.02 of 0.  (Hint: There is a solution with F1 = 0 and F3 = 0 and F4 = -.09 so you only need to fiddle with F2 to find an appropriate number.)

The goal of this exercise is not to design a control system but to illustrate that eigenvalues and eigenvectors are associated with basic behavior (called resonant modes) of the airplane. This is true not just for mechanical objects but for anything which changes with time. Learning one example well, such as this airplane situation, can give intuition for a large number of areas.

§4.5 Feedback

Introduction to MATLAB
We're almost done, but before we finish we'd like to get some (optional) feedback from you so that we can improve this lab.  If you have any comments, write them (by hand if you like) on a separate page at the end of your assignment, so that we can tear it off before we give back your assignment.

Feedback (Optional)
Do you have any suggestions to help improve this lab?  This includes errors or potential errors, typos, unclear statements or questions, or anything else that comes to mind.  The more specific you can be, the easier it will be for us to improve the lab.

§4.6 Conclusion

    Systems of differential equations consitute the math model central to in many technological and scientific applications.  In practice, models requiring many differential equations are much more common than models using only one. Dozens of variables are common.   In this lab we saw how matrices and a little linear algebra can give us powerful tools for working with linear systems, even very large ones.  

Later in class you will study Laplace transforms. For linear systems they combine very well with the linear algebra techniques we have seen here to produce some of the main design techniques used in engineering.

Last Modified:

11/02/2009