How to plot isochrons in 3d
3 views (last 30 days)
Show older comments
I have code for plotting isochrons for a system of ODE's in 2 dimensions, but I need to plot isochrons for a system of 3 ODEs in 3 dimensions. I have already made some adjustments to the plot commands and such, but it is not running the way I intended. Any help or advice would be greatly appreciated !
[t,lc,iso]=isochron
function [t,lc,iso]=isochron()
alpha1 = 0.1;alpha2 = 3;alpha3 = 3;beta1 = 3;beta2 = 1;beta3 = 1;K1 = 0.5;K2 = 0.5;K3 = 0.5;n1 = 8;n2 = 8;n3 = 8;
system = @(t,Y)[alpha1 - beta1*Y(1)*(Y(3)^n1/(K1^n1+Y(3)^n1));
alpha2*(1-Y(2))*(Y(1)^n2/(K2^n2+Y(1)^n2))-beta2*Y(2);
alpha3*(1-Y(3))*(Y(2)^n3/(K3^n3+Y(2)^n3))-beta3*Y(3)];
T=1.395*pi; % period of the oscillator
%x0 = [0.226623;0.174165;0.243654];
x0=[0.455676;0.555023;0.504329]; % initial condition (must be on the limit cycle)
num_chron=5; % number of isochrons plotted
phases=0:T/num_chron:T;
tau=T/2000; % time step of integration
m=200; % spatial grid
k=0; % the number of skipped cycles
[t,lc] = ode45(system,0:tau:T,x0); % call the DE
dx=(max(lc)-min(lc))'/m; % spatial resolution
center = (max(lc)+min(lc))'/2; % center of the limit cycle
iso=[x0-m^0.5*dx, x0+m^0.5*dx]; % isochron?s initial segment
for t=0:-tau:-(k+1)*T % backward integration
for i=1:size(iso,2)
iso(:,i)=iso(:,i)-tau*feval(system,t,iso(:,i)); % move one step
end
i=1;
while i<=size(iso,2) % remove infinite solutions
if any(abs(iso(:,i)-center)>1.2*m*dx) % check boundaries
iso = [iso(:,1:i-1),iso(:,i+1:end)]; % remove
else i=i+1;
end
end
i=1;
while i<=size(iso,2)-1
d=sqrt(sum(((iso(:,i)-iso(:,i+1))./dx).^2)); % normalized distance
if d > 2 % add a point in the middle
iso = [iso(:,1:i), (iso(:,i)+iso(:,i+1))/2 ,iso(:,i+1:end)];
end
if d < 1 % remove the point
iso = [iso(:,1:i), iso(:,1:i)];
else i=i+1;
end
end
if (mod(-t,T)<=tau/2) && (-t<k*T+tau)
cla;
plot3(lc(:,1),lc(:,2),lc(:,3),'r'); % plot the limit cycle
hold on;
end
if min(abs(mod(-t,T)-phases))<tau/2 % plot the isochrons
plot3(iso(1,:),iso(2,:),iso(3,:),'k-');
disp(i<=size(iso,2)-1)
drawnow;
end
end
xlim(0:1)
ylim(0:1)
zlim(0:1)
hold off;
end
0 Comments
Answers (1)
Arun
on 17 Oct 2023
Hey Braden,
I understand that you want a 3-D plot of “isochrons” for a 3-ODEs and have encountered difficulties in obtaining the desired outcomes.
Here is a way to obtain 3-D plot:
if (mod(-t,T)<=tau/2) && (-t<k*T+tau)
cla;
scatter3(lc(:,1),lc(:,2),lc(:,3),'r'); % plot the limit cycle
hold on;
end
if min(abs(mod(-t,T)-phases))<tau/2 % plot the isochrons
scatter3(iso(1,:),iso(2,:),iso(3,:),'k');
disp(i<=size(iso,2)-1)
drawnow;
end
Output plots for your reference:
Hope this helps.
0 Comments
See Also
Categories
Find more on Digital Filter Analysis 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!