Lab 8 -- Lorentz Systems
Lab 8 (click here for pdf) concerns solving the system
x' = s(y - x)
y' = rx - y - xz
z' = -bz + xy
Outside of lab: Just remember to do this part. It shouldn't be any problem.
- Follow the steps to create a function file for the Lorentz system corresponding to s = 10, r = 28, and b = 8/3. Enter the lines:
----
function up = lor(t,u)
up(1,1)=-10*u(1) + 10*u(2);
up(2,1)=28*u(1)-u(2)-u(1)*u(3);
up(3,1)= u(1)*u(2)-8/3*u(3);
----
Where "up," a column vector, corresponds to the vector u', and u is your input vector (u = [u(1), u(2), u(3)] = [x, y, z], which makes sense of the system). Use the command
[t, u]=ode45('lor', [0, 20], [0;1;0]); plot3(u(:,1),u(:,2),u(:,3));
and witness the butterfly shape. We're going to find formulas for the planes that the "wings" of the butterfly shape approaches.- To find the behavior around the hole in the "wings," we want to linearize the plane around the center of that hole, which is the equilibrium point (6*2^0.5,6*2^0.5, 27). We want to figure out what the system behaves like around this point.
We need the linearized system to be centered at the equilibrium point to draw conclusions about the behavior around that equilibrium point. There are two steps: (1) Shift the system and (2) linearize the system.
Step 1: To shift the system, we let [u, v, w] = [x, y, z] - (the equilibrium point). So [u, v, w] = [x-6*sqrt(2), y-6*sqrt(2); z - 27]. Thus [u, v, w] = [0, 0, 0] corresponds to [x, y, z] = (the equilibrium point).
You can now write u' = ( x-6*sqrt(2) )' = 10(y - x) = 10( (v+6sqrt(2)) - (u+6*sqrt(2) ) = 10(v - u). Now find v' and w' in terms of u, v, and w.
Step 2: Linearize the system. Now you can write it in the form X' = AX, where X = [u, v, w].- Use Matlab to get the eigenvectors and eigenvalues.for the matrix
A = [-10 10 0; 1 -1 -6*2^.5; 6*2^.5 6*2^.5 -8/3 ];
Use the commands given to get the eigenvalues. Get the real and imaginary components of the complex eigenvectors as in the instructions.- After you plot the system, the tricky part will be to prove that the orbits lie on the same plane. Let's start by giving the general solution to a system of the form X' = AX where A is 3x3.
x(t) = c1*v1*e^(lambda_1*t) + c2*v2*e^(lambda_2*t) + c3*v3*e^(lambda_3*t) is the general solution. But what does this have to do with solving the system x' = Ax in terms of R and I?Two eigenvectors of A are complex. You will find that the real component of one complex eigenvector, R, is a linear combination of the other complex eigenvector's real component. The complex component of that eigenvector, I, is also linearly dependent on the complex component of the other eigenvector.
Well, two of the eigenvectors, v2 and v3 are complex. Remember that the real and complex parts of the eigenvectors are linearly dependent. You will notice that v2 = R + i*I and v3 = R - i*I.
Now let m(t) = c2*e^(lambda_2*t) and let n(t) = c3*e^(lambda_3*t).
x(t) = c1*v1*e^(lambda_1*t) + m(t)*(R+i*I) + n(t)(R-i*I)
= c1*v1*e^(lamdba_1*t) + (m(t)+n(t))*R + i*(m(t)-n(t))*I.Let f(t) = m(t) + n(t) and let g(t) = i*(m(t) - n(t)).
Then x(t) = c1*v1*e^(lambda_1*t) + f(t)*R + g(t)*I.
And this ought to put us all at ease because now
(1) R, I, and v1 are linearly independent and
(2) We solve the initial conditions
x(0) = R
R = x(0) = c1*v1*e^(0) + f(0)*R + g(0)*I
= c1*v1 + f(0)*R + g(0)*I
==> c1 = 0, f(0) = 1, g(0) = 0
(Why is this implied?)
x(0) = I
I = x(0) = c1*v1*e^(0) + f(0)*R + g(0)*I
==> c1 = f(0) = 0, g(0) = 1.
(Why is this implied?)Thus the solutions to X' = AX with x(0) = R and x(0) = I fall in a single plane, RxI (why?). That is to say, the span of two linearly independent eigenvectors R and I is the desired plane.
The plane in question is the span of R and I. To find an equation of the form au + bv + cw = 0, consider that RxI is perpendicular to the plane in question. Now note that the plane (call it P) that is perpendicular to RxI, is described by (RxI) dot P = 0.- Let W = the real eigenvector from P (P comes from #3). Hint for proving that the orbit is a line: relate your answer to word "span." Find the limit of your solution as t -> infinity.
- Clear plot and plot the solution with initial data (R+seed/10*W)). Explain what you think to be the behavior of the orbits with respect to the plane RxI.
- Enter "cla" and "hold off." Plot the phase portrait for the original (non-linearized) Lorentz system with the given intervals.
- When thinking about this one, it might help you to notice (1) the equilibrium point to the original system, and (2) a correspondence between between T + rX + sY and the linearized system from #2. Note that we analyzed the behavior of a linearized version of this system with the equilibrium point shifted to (0,0,0) in #4. Here our equilibrium point is not (0,0,0) but T.
- Take the equation [u,v,w] dot (RxI) that you used earlier, but re-write u, v, and w in terms of x, y, and z. (Remember what you did in #2.)
- To really do this one to full length, you would first shift your original Lorentz system to be centered around the equilibrium point at the center of the other plane (-6*sqrt(2), -6*sqrt(2), 27). You would then linearize the shifted system, and find the R and I for this plane. (It does not turn out that the R and I are the same as for the original system.) The plane that the other "wing" is attracted to is in the form T + rX + sY. You would use the same process as you used in number 9 (but with your adjusted values of u, v, and w).
- Plot your original Lorentz system. Try view(a,b) a few times and print out two plots that reveal different aspects of the system. I.e., you're in good shape if your two printouts don't look exactly the same.
- Follow the instructions for this one. We'll only ask that you print out at least one plot for part (a) and one plot for part (b). These should be 2d plots.
That should be all for this lab (good riddance).
|
e-mail: renn@purdue.edu |