3d phase portrait

Hi, I want to plot a 3-d phase portrait for a set of 3 ODEs, i have used the following code and i get a basic plot, i was wondering how to add direction arrows and a mesh grid and why i only get single spirals for the trajectories. Here is the code i have used:

function f = cluster(t,y)
  %BD
  a=1;
  b=1.2;
  %equilibrium values
  %c1 equilbrium value

    f = zeros(size(y));
    f(1) = -50*a*y(1)-b*y(1)+15*a*y(1)*y(2)+20*a*y(2)*y(3)+y(2)*b+9*a*y(2)^2+6*a*y(1)^2-60*a*y(2)-80*a*y(3)+24*a*y(2)*y(3)+16*a*y(3)^2;
    f(2) = 10*a*y(1) - a*y(1)*y(2) -4*a*y(1)*y(3) -2*a*y(1)^2 -b*y(2) +3*a*y(2)^2 -10*a*y(2) +4*a*y(2)*y(3) +b*y(3);
    f(3) = -2*a*y(1)*y(2) - 3*a*y(2)^2 -4*a*y(2)*y(3) +10*a*y(2)-b*y(3);

Plotting code:

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

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

I have computed the corresponding eigen values and vector points for specific equilbrium points in a separate file not sure if that will help?

[Akshay Khadse]

Can you elaborate on what is the "f" in your code above? According to me, for a phase portrait, "f" should be the gradients. However, you are plotting the solution of the differential equations, hence the single spirals.

[Gulmira Tussupbekova]

Can you please send a code for this problem

Answer 1

[Akshay Khadse]

You can get 3D Phase Portraits by plotting the gradients against the co-ordinates using the " meshgrid ", and " quiver3 " functions.

" meshgrid " is used to generate the 2D or 3D grids and " quiver " or " quiver3 " is used to place arrows at these co-ordinates.

Creating meshgrid:

[x1,y1,z1] = meshgrid(-2:0.2:2,-2:0.2:2,-2:0.2:2);
u = zeros(size(x1));
v = zeros(size(y1));
w = zeros(size(z1));

Calculating gradients:

t=0;
for i = 1:numel(x1)
    Yprime = cluster(t,[x1(i); y1(i); z1(i)]);
    u(i) = Yprime(1);
    v(i) = Yprime(2);
    w(i) = Yprime(3);
end

Plotting:

quiver3(x1,y1,z1,u,v,w); figure(gcf)

Please refer to the documentation for "quiver3" here for examples.

[Aishah Malek]

This makes sense , thankyou for you help