How can i draw the phase portrait in my code for 3 non-linear coupled ODEs in 3D plane

function SC1C2

clear

clc

clf

close all

%==== parameters =====

k1 = 1.5;

k_1 = 0.5;

k2 = 1;

k3 = 2;

k_3 = 1;

k4 = 3;

%======Initial conditions =====

initialx = 4;

initialy = 0;

initialz = 0;

runtime = 10; % Simulation run time

iterations = 1; % Sets iteration no. to 1

pausetime = 0.01; % Pause time between animation

az=-37.5; % Sets azemouth angle

el= 30; % Sets elevation angle.

deq1=@(t,x) [k_1*x(2)-k1*x(1)-k3*x(1)*x(2)-k_3*x(3);

k1*x(1)-(k_1+k2)*x(2)-k3*x(1)*x(2)+(k4+k_3)*x(3);

k3*x(1)*x(2)-(k4+k_3)*x(3)];

[t,sol] = ode45(deq1,[0 runtime],[initialx initialy initialz]);

arraysize = size(t);

for i = 1 : max(arraysize)

plot3(sol(iterations,1),sol(iterations,2),sol(iterations,3),'-o','markersize',6,'MarkerFaceColor','m')

hold on

grid on

xlabel('S','fontsize',12)

ylabel('C1','fontsize',12)

zlabel('C2','fontsize',12)

iterations = iterations + 1;

pause(pausetime)

axis([min(sol(:,1)) max(sol(:,1)) min(sol(:,2)) max(sol(:,2)) min(sol(:,3)) max(sol(:,3))])

view(az,el)

az = az+1;

end

end

Answer 1

[William Rose]

@Akhtar Jan,

Matlabs quiver3() is very nice for this.

"quiver3(X,Y,Z,U,V,W) plots arrows with directional components U, V, and W at the Cartesian coordinates specified by X, Y, and Z. For example, the first arrow originates from the point X(1), Y(1), and Z(1), extends in the direction of the x-axis according to U(1), extends in the direction of the y-axis according to V(1), and extends in the direction of the z-axis according to W(1). By default, the quiver3 function scales the arrow lengths so that they do not overlap."

[Akhtar Jan]

Thanks for reply william actually i tried it but couldnt draw its in 3d i have 3 variables

[Akhtar Jan]

Can u pls edit in the code?

[William Rose]

@Akhtar Jan,

This shows code you can add to your function, after the for loop which is at the end of your function. This code plots an arrow at each simulation solution point, except the last point. Each arrow points in the direction of travel. In other words, each arrow points toward the next solution point.

function SC1C2
%code here...
for i = 1 : max(arraysize)
    %more code here...
end
Xarrow=sol(1:end-1,1);      %x locations for arrows
Yarrow=sol(1:end-1,2);      %y locations for arrows
Zarrow=sol(1:end-1,3);      %z locations for arrows
Uarrow=diff(sol(:,1));      %x length of arrows
Varrow=diff(sol(:,2));      %y length of arrows
Warrow=diff(sol(:,3));      %z length of arrows
quiver3(Xarrow,Yarrow,Zarrow,Uarrow,Varrow,Warrow); %plot arrows
end     %end of function SC1C2

Try it. The arrows are hard to see. The default arrow color on my machine is an orange color which does not show up well. You can make the arrows black as follows:

quiver3(Xarrow,Yarrow,Zarrow,Uarrow,Varrow,Warrow,'k'); %plot black arrows

The arrows are still difficult to see, because each arrow bumps into the next point.

[Akhtar Jan]

i have run it thanks let me ask one question Is there any way to make the arrows size large? also can we draw the arrows along the direction with the plot i mean in the whole x,y,z plane?

Regards

[William Rose]

@Akhtar Jan,

"Is there any way to make the arrows size large?"

Yes. See options AutoScale and AutoScaleFactor, described in the help for quiver3. You may turn off AutoScale. Then the arrows have the legths specified in vectors U,V,W, and you could pass scaled versions of them:

sf=1.5;
quiver3(x, y, z, sf*u, sf*v, sf*w, 'AutoScale','off');

or you could leave Autoscale on, but change the AutoScaleFactor to a number larger than the default, which is 0.9. See the help.

"can we draw the arrows along the direction with the plot i mean in the whole x,y,z plane?"

I assume that you are asking if it is possible to draw arrows throughout the x,y,z volume. Yes, it is possible. In what direction should each arrow point? In the direction of the derivatives at that point. Your function deq1() computes the three components of the derivative. Fortunately, it does not depend on time. Therefore the direction of the derivatives is time-invariant. In the code below, I compute and plot arrows at 8x8x8 locations.

function SC1C2
%code here...
for i = 1 : max(arraysize)
    %more code here...
end
N=8;        %number of arrows along each dimension
%Next: Compute N evenly spaced points along each axis.
xvals=linspace(min(sol(:,1)),max(sol(:,1)),N);
yvals=linspace(min(sol(:,2)),max(sol(:,2)),N);
zvals=linspace(min(sol(:,3)),max(sol(:,3)),N);
%Next: Make vectors of length N^3 that contain the x,y,z locations of all the arrows.
[X,Y,Z]=meshgrid(xvals,yvals,zvals);  %creates NxNxN arrays
%Next: Make X,Y,Z row vectors, since quiver3() won't accept NxNxN array.
X=reshape(X,1,N^3);
Y=reshape(Y,1,N^3);
Z=reshape(Z,1,N^3);
U=zeros(1,N^3); V=U; W=U;  %pre-allocate arrays U,V,W
for i=1:N^3
    %Next: Compute the direction components [U,V,W] for each arrow.
    deriv=deq1(0,[X(i),Y(i),Z(i)]);
    U(i)=deriv(1); V(i)=deriv(2); W(i)=deriv(3);
end
quiver3(X,Y,Z,U,V,W,'k'); %plot black arrows

end  %end function SC1C2

Try it.

[Akhtar Jan]

Thank you sir....its run.👍👍👍❤❤❤

[Akhtar Jan]

now can i plot the first graph among these arrows i mean the phase portarit of my odes?

[William Rose]

@Akhtar Jan,

"can i plot the first graph among these arrows i mean the phase portarit of my odes?"

We can also include the arrows along the path, by adding back the first code fragment. Let's use blue arrows to make them more visible. The lengths of the black and blue arrows may be scaled differetly, due to the autoscaling. Code attached.

[Akhtar Jan]

thanks alot sir it works brilliant you are the greatest sir i love you..........

may God give you long and healthy life.ameen

[William Rose]

@Akhtar Jan, thank you very much for you kind comments. Good luck with your work!