Problem regarding a periodic solution not being saved properly

Hello! I am working on a problem where I have to solve an ODE for time of 10 years and then check when and where that solution becomes periodic. I want to store the value of the periodic solution and the interval when it becomes periodic. This is what I tried so far:

clear all;

p=10;% Number of year for which simulation runs
y0=20000; % Initial condition
tspan=0:1:10*365;
epsilon=1e-10; % Tolerance for periodicty condition which is defined further
t=0:3650; % This t will be used in the following parameter m files
%Calling some parameters
mugen;
Kgen;
d1gen;
delta1gen;

[X Y] = ode15s(@(t,y)para1(t,y,mug,Kg,d1g,delta1g),tspan,y0);

if (X>=365) & (abs(Y(X)-Y(:,X-365))<epsilon) % Check for every
                                              %time value if the solution 
                                              %is periodic by using the definition of periodicity.
   YcSteady(X)=Y(:,X); % Store the periodic the solution found
   XcSteady(X,:)=[X-365  X]; % Store time interval when periodic solution occurs
end

Now, the problem is that it does not store YcSteady which is the time at which periodicity condition is satisfied and XcSteady which is the time interval between which the condition is satisfied!

Thanks!

[Star Strider]

It would be easier to answer if we could run your code. We have no idea what X and Y are numerically, so we cannot evaluate your logic. We at least seem to be missing para1. Attaching it would help.

[Rose]

Attached is the file para1.m. Is there anything else that i should attach?

[Rose]

Attached!

[Rose]

Attached are the parameter files that I am calling in the code. Thanks for your time.

Answer 1

[Star Strider]

I’m still not certain exactly what you want to do, but I hope to have made it easier for you to do it. I’m only posting part of the main code (the part I added to). I didn’t alter anything else.

The added code subtracts the mean from Y to create zero-crossings (that are actually mean-crossings), then detects them and reports their indices as the Xzx variable. It takes every other zero-crossing as the definition of a period, and reports their X-coordinates in the Pds variable. (The Prs variable calculates the mean and standard deviation on Pds.) The plots in figure(1) are X,Y, the mean, and the mean-crossings. I did not change the values of your X or Y data.

The Y value does not cross the mean line until the onset of the first complete period at Pds(1), X = 254, so the onset of periodicity would seem to me to be Pds(1)), and the duration of periodicity to Pds(end). I will defer to you as to how you want to define periodicity. You can probably define XcSteady and YcSteady from the data I calculated. If you want to define the onset of periodicity differently with my code, the easiest way is to change mean(Y) to Thld, where Thld is the scalar value of your choice. (It’s best not to change anything else in my code.) I don’t understand XcSteady and YcSteady, so I did not change any part of your code. I simply added mine. Experiment with it to understand what it does.

The code:

[X Y] = ode15s(@(t,y)para1(t,y,mug,Kg,d1g,delta1g),tspan,y0);

% Find mean crossings: 
Ymz = Y - mean(Y);                          % Subtract mean creating zero-crossings 
Ymx = Ymz .* circshift(Ymz, [-1 0]);        % Negatives are zero-crossings
Xzx = find(Ymx <= 0);                       % Indices of zero-crossings
Pds = X(Xzx(1:2:end-1));                    % Define periods as values of X
Prs = [mean(diff(Pds))  std(diff(Pds))];    % Statistics on the periods

if (X>=365) & (abs(Y(X)-Y(:,X-365))<epsilon) % Check for every
                                              %time value if the solution 
                                              %is periodic by using the definition of periodicity.
   YcSteady(X)=Y(:,X); % Store the periodic the solution found
   XcSteady(X,:)=[X-365  X]; % Store time interval when periodic solution occurs
end

figure(1)
plot(X, Y)
hold on
plot(X(Xzx), Y(Xzx), '*r')
plot(X, ones(size(X))*mean(Y), '-g')
hold off
grid

[Rose]

Thank you very much for such a neat and detailed reply. Though it was not exactly what I needed, but this was entirely new thing that I learnt. The plot was so good and informative. I really appreciate that you spent your precious time for writing solution to this problem. Thanks once again..!

[Star Strider]

My pleasure!