Problem regarding a periodic solution not being saved properly

2 次查看（过去 30 天）
Rose
Rose 2014-7-5
评论： Star Strider 2014-7-8
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!
  4 个评论
显示 2更早的评论
Rose
Rose 2014-7-5
Attached are the parameter files that I am calling in the code. Thanks for your time.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2014-7-5
编辑：Star Strider 2014-7-5
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
  2 个评论
Rose
Rose 2014-7-8
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..!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by