Poincare map of Lorenz equations (homework)
23 views (last 30 days)
Show older comments
So I am working on a homework assignment and the first part was to plot the lorenz equations on a 3d map. Which I believe I did. Then I needed to modify the code for every time the trajectory crosses the plane z = 27, plot the corresponding point in the xy-plane. And do this for the time interval [30, 70]. Thats where my code gets messed up. This is what I have:
t(1)=0;
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10;
rho=28;
beta=(8/3);
h=0.01;
t=0:h:70;
f=@(t,x,y,z) sigma*(y-x);
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1)
k1=f(t(i),x(i),y(i),z(i));
l1=g(t(i),x(i),y(i),z(i));
m1=p(t(i),x(i),y(i),z(i));
k2=f(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
l2=g(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
m2=p(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
k3=f(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
l3=g(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
m3=p(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
k4=f(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
l4=g(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
m4=p(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
x(i+1) = x(i) + h*(k1 +2*k2 +2*k3 +k4)/6;
y(i+1) = y(i) + h*(l1 +2*l2 +2*l3 +l4)/6;
z(i+1) = z(i) + h*(m1+2*m2 +2*m3 +m4)/6;
end
figure(1)
plot3(x,y,z)
for i=30:70 ****************This is where the code gets messed up**************
for z=27
x(i)=x;
y(i)=y;
end
figure(2)
plot(x(i),y(i))
end
My thought process was to use a for loop first with the time interval condition then with the condition that z=27. Then take x(i) which was solved in the runge kutta code and as the value of i gets plugged in when the z value is 27, it would give an x coordinate and y coordinate for when the plane is crossed at z=27. But doing this loop messes up the code to where I cant even plot figure 1 anymore. I feel I am using the for loop wrong, but I thought doing something like this would work. Any hints or documentations to read to help me get on the right track would be greatly appreciated.
0 Comments
Accepted Answer
Image Analyst
on 2 Nov 2021
Not sure what this means or what your intent is:
x(i)=x;
y(i)=y;
x and y are 7001 long vectors. And you're trying to stuff the whole vector into just one single one of their elements.
If you want to find indexes of x, y, and z where z is 27, then create a mask for where z is in that tight range:
mask = z > 26.9 & z < 27;
% Then extract the x and y values:
x27 = x(mask);
y27 = y(mask);
% Then plot those
plot(x27, y27, 'b-');
0 Comments
More Answers (0)
See Also
Categories
Find more on 2-D and 3-D Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!