18.03 Problem Set 6

Due by 1:00 P.M., Monday, April 12, 1999, in the boxes at 2-106, next to the Undergraduate Mathematics Office.

Syllabus

23. (F 2 Apr) Complex or repeated eigenvalues: Notes LS 1.3, 1.4.

24. (M 5 Apr) Dynamics of linear autonomous systems: Notes GS 1–5, EP 7.3 (522– 530), Portrait Gallery Handout.

25. (W 7 Apr) Dynamics of nonlinear autonomous systems; linearization: Notes GS 6–7.

26. (F 9 Apr) Example: the nonlinear pendulum: EP 7.5 (561).

27. (M 12 Apr) Initial Value problems: Notes LS 2.1–2.3; PS6 due.

28. (W 14 Apr) Matrix exponentials: EP 5.7 (441–443).

29. (F 16 Apr) Spinning books—a review: Handout on Euler’s equations.

30. (W 21 Apr) Hour Exam III

Part I.

23. (F 2 Apr) EP 5.4: 1, 5, 13, 17.

24. (M 5 Apr) Notes NL: 8.

25. (W 7 Apr) Notes NL: 13, 14.

26. (F 9 Apr) Use pplane5 to generate a phase portrait of the damped nonlinear pendulum θ̈ + θ̇ + 2 sin(θ) = 0, 0 ≤ θ ≤ 12, −2 ≤ θ̇ ≤ 2. Hand in a printout.

Part II.

23. (F 2 Apr) This problem will analyze the family of first order homogeneous linear differential equations

\{ \begin{matrix}ẋ & = & y \\ ẏ & = & −ax − 2y. \end{matrix}

where a is a constant. This is the system associated to the second order equation ẍ + 2ẋ + ax = 0.

(a) What is the characteristic polynomial of the matrix A of coefficients of this system? Compute the eigenvalues λ_{1}, λ_{2}. For each eigenvalue, find a nonzero eigenvector. For what value(s) of a is there a single (repeated) eigenvalue? Write down the matrix A in that case. Does it have two independent eigenvectors? For what values of a are the eigenvalues real and distinct? For what values of a are the eigenvalues nonreal complex numbers? Note that they indeed complex conjugates of one another in this last case. For what value(s) of a is 0 an eigenvalue? What is the other eigenvalue in that case?

(b) Write down the general real solution to the system ˙ u⃗ = Au⃗. You will have to separate into three cases, according to whether the roots are real and distinct, equal, or nonreal. The hardest case is when there is repeated eigenvalue. Here is a hint for that case: if α⃗_{1} is a nonzero eigenvector, so that we have one solution of the form e^{λ_{1}t}α⃗_{1}, find α⃗_{2} so that e^{λ_{1}t}(tα⃗_{1} + α⃗_{2}) is another solution.

24. (M 5 Apr) This problem will use the Matlab tool pplane5. to understand the qualitative behavior of the solutions you found in 23. So enter Matlab and type pplane5. In the PPLANE5 Setup window, remove the default equations and enter y for x' and -a*x-2*y for dy/dt. In the first window following “Parameters:” enter a and then -8. Your plots will look better if you reset the maximal value of x to 2 and the minimal value of y to -2. Do this and type proceed. A new window labelled PPLANE5 Display will appear, with the direction field nicely displayed. If you position the cursor in the window and tap the left mouse key a trajectory will form. Do this a dozen times, for various starting positions, to generate a phase portrait. Be sure to begin several trajectories near the origin.

(a) Rather than printing this out, carefully copy this phase portrait onto your answer sheet. Then return to the PPLANE5 Setup window and change the value of a to -1. Repeat the process, and then do it again, till you have phase portraits for a = −8, −1, 0, .5, 1, 2, and 10. Label each with the corresponding value of a. Line them up so you can see the progression of shapes as a increases from −8 to 10. Label each: is it a saddle, a proper node, an improper node, a star node, a stable center, a spiral, or none of the above? Indicate the direction of travel along the trajectories. Is the equation asymptotically stable, neutrally stable, or unstable?

(b) For a < 1 (but a ≠ 0) some of the trajectories should be almost straight rays

eminating from the origin. The picture hints that there are solutions u⃗(t) = [ \begin{matrix}x(t) \\ y(t) \end{matrix} ]

which parametrize such rays. Is this the case? Find them analytically. What are the slopes (in terms of a, of course)? Why did we actually know in advance that in this case some solutions would lie along rays from the origin?

(c) When a = 0 the trajectories should all look like parallel straight rays. Are they? What is the slope? In this case there are some solutions missing from the work pplane5 did for you. These are constant solutions, whose trajectories are just points. You are unlikely to have chosen them as starting points. By looking at the analytic solution, find these constant solutions. What feature of the matrix lead to the existence of so many constant solutions? What constant solutions exist for any first order homogeneous linear system?

(d) When a = 1 your picture should indicate a single pair of opposed rays as solutions. What is the slope of those solutions? Draw a new picture showing the trajectories of the two independent solutions you found for this case above. Then show that with the exception of the solutions along the rays and the constant solution u⃗ = 0, every solution crosses both the x-axis and the y-axis exactly once. (This is hard to see on the picture, for some trajectories.)

Your picture for a = 2 should look quite a bit like the picture for a = 1. This is actually deceptive, though. When a = 1, you saw that x (and y) change sign at most once, while when a = 2 the system is underdamped and the solutions change sign infinitely often. You don’t see this in ppland5 because the spirals are so small or so large. (You can go back to this value of a and zoom in towards the origin using Edit followed by Zoom in, if you like. There is a fractal quality to these pictures; they look the same under magnification as they did at first.) As a increases they grow, and you can see more of them in the a = 10 picture.

25. (W 7 Apr) (This problem is not directly relevant to the lecture of that day.) We’ll use Matlab to get a feel for the solutions of the nonlinear jump rope y′′+y\sqrt{1 + y^{′2}} = 0. Matlab numerical ODE routines don’t handle second order equations. This is because any second order equation can be expressed as a system of first order equations. To do this, we denote y′ by a new variable name, say z(x). Then z′ = y′′ = −y\sqrt{1 + y^{′2}} = −y\sqrt{1 + z^{2}}: we have a nonlinear first order system

\{ \begin{matrix}y′ & = & z \\ z′ & = & −y\sqrt{1 + z^{2}}. \end{matrix}  (1)

The Matlab function ode45 can handle vector-valued ODEs just as easily as ordinary ones. I need to say a word about how Matlab treats vectors and matrices. You can creat a row-vector by typing for example r=[2,2] or equivalently [2 2]. Try it. Then type c=[2;2]: you’ll get a column vector. A general matrix can be created by combining these: mat=[1 2;3 4]. The nth entry in a row-vector or a column vector is obtained by for example r(2) or c(1). The entry in the first row and second column of mat is mat(1,2). If you type merely mat(2) you’ll get the first entry in the second row, strangely enough. If you want to extract the entire second row you need to type mat(2,:). The first column is mat(:,1). Try it.

Matrix multiplication is denoted just as multiplication of numbers, with a star. Try r*c, c*r, and mat*mat. Do the answers make sense? What happens if you try c*c or r*r? Powers are formed using ^. Try mat^10 for example.

Determinants are computed using det. What is det(mat)? How about det(mat^10)? Does this accurately reflect the wonderful fact that det(AB) = det(A) det(B) (if A and B are square matrices of the same size)?

Now create a new function in a file called jrope.m, in the directory from which you launched Matlab, containing the following:

function yprime=jrope(x,y)
yprime=[y(2);-y(1)*sqrt(1+y(2)*y(2))];

In this function y is a column matrix with two entries: the first gives values of the function y in (1), and the second gives values of the function z in (1). yprime is also a column vector with two entries, giving y′ and z′. ode45 will use this function to produce a column-vector of x-values, and a matrix with two columns: the first is the list of yvalues corresponding to the given x-values, and the second is the list of z-values. This is confusing: the rows in the second output matrix correspond to the columns in the function.

So now type

[x,y]=ode45('jrope',[0,10],[0,1])

and watch what comes out. You can plot the y-values, which are the ones we care about, by writing

plot(x,y(:,1)).

Let’s add a grid by typing grid. This plot should look pretty sinusoidal; it’s pretty well approximated by the linear model. Now let’s increase the initial slope to 4:

[x,y]=ode45('jrope',[0,10],[0,4])

Before plotting it, we’d better type hold on so the original graph isn’t replaced. Then type plot(x,y(:,1),'r'). Do this (using 'g', 'k', 'm') for y′(0) = 16, 64, and 256. Print out the result (black and white is fine) and hand in the result as (a).

The jump rope handout indicates that the amplitude of y is given by the formula y_{max} =

\sqrt{2 (\sqrt{1 + y′(0)^{2}} − 1)}. For y′(0) = 256 this gives approximately 22.5832660660558,

which should be pretty close to what your graphs show.

Notice how unlike the linear case this is. Solutions to y′′ + y = 0 with y(0) = 0 all have the form y′(0) sin(x). They all have the same period, 2π, and the amplitude equals |y′(0)|. On the other hand, the solutions to the nonlinear equation y′′ + y\sqrt{1 + y^{′2}} = 0 with y(0) = 0 have period shortening with y′(0) (but beginning close to 2π when y′(0) is small) and amplitude growing roughly like the square root of |y′(0)|. (Right?—according to the graphs, when you multiply y′(0) by 4 the amplitude roughly doubles.)

Notice also that as y′(0) grows the solutions look less and less sinusoidal: the tops are progressively flatter. I guess this is good; it gives you more room to jump (although it doesn’t entirely make up for the shortening period).

We also may care about the tension along the rope. As explained in the handout, this is given by \sqrt{1 + y^{′2}} (times a constant T which for convenience I’ll take to be 1). This is close to |y′| when |y′| is large. We’ll plot it. Create a matrix y with first column a sequence of values of y and second column a sequence of values of z = y′, for y′(0) = 256,

from x = 0 to x = 3, by typing [x,y]=ode45('jrope',[0,3],[0,256]). Plot this by

killing Figure No. 1 and then typing plot(x,y(:,1)). Prepare to plot the tension function as well by typing hold on. Create a new vector whose entries are corresponding tension values: tension=sqrt(1+x(:,1)*x(:,1)). You should get an error message. You have tried to multiply two column vectors together, and you can’t do that. You want to multiply their entries, one by one, and create a new column vector. You can do that by typing .* instead of *: so try tension=sqrt(1+x(:,1).*x(:,1)) and then plot(x,tension,'r'). How does it look? You can see it better by typing grid. The tension should always be positive, bigger than or equal to 1 in fact, and largest where the slope of the jump rope is greatest. Print out the result and turn it in as (b).