Homework 1 Matlab Help

To make 3D phase portraits in matlab we use the generic ODE solvers, particularly the ode45 routine. First create a file called lorenz.m which contains the following code to evaluate the vector field at a given time t and position y.

function f = lorenz(t,y)

sigma = 10;
r = 28;
b = 8/3;

f = zeros(size(y));
f(1) = sigma * (y(2) - y(1));
f(2) = r*y(1) - y(2) - y(1)*y(3);
f(3) = y(1)*y(2) - b*y(3);

Next we can generate a trajectory of times t and phase-space points y by running the following function inside matlab (note that it will take some minutes for this command to finish).

[t,y] = ode45(@lorenz,[0:0.01:100],[1 2 3]);

Here the arguments say to use the lorenz.m file for the vector field, to use the time seqence [0:0.01:100] (start at t = 0, end at t = 100 and use a timestep of 0.01), and to take the initial condition to be the point [1 2 3].

The trajectory generated from this initial condition is stored in the vector t of times and the matrix y of phase-space positions. The vectors y(:,1), y(:,2), y(:,3) give the (x,y,z) vectors of the trajectory. We can plot different views of the trajectory with commands such as

plot(t,y(:,3));                 % plot of z(t) versus time
plot(y(:,1),y(:,3));            % plot of z versus x
plot3(y(:,1),y(:,2),y(:,3));    % 3D plot of trajectory

To see the trajectory being generated in real-time, copy the file linephas3.m and then use the commands:

options = odeset('OutputFcn',@linephas3);
[t,y] = ode45(@lorenz,[0:0.01:100],[1 2 3],options);

To compute the difference between two nearby trajectories, you can use commands such as the following.

[t,ya] = ode45(@lorenz,[0:0.01:100],[1 2 3]);
[t,yb] = ode45(@lorenz,[0:0.01:100],[1.0000000001 2 3]);
delta = abs(ya(:,1) - yb(:,1));
plot(t,log(delta));

[t,ya] = ode45(@lorenz,[0:0.01:100],[1 2 3]);
[t,yb] = ode45(@lorenz,[0:0.01:100],[1.0000000001 2 3]);
plot(t,ya(:,1));
hold on
plot(t,yb(:,1));