Problem regarding a periodic solution not being saved properly
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Rose
am 5 Jul. 2014
Kommentiert: Star Strider
am 8 Jul. 2014
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 Kommentare
Akzeptierte Antwort
Star Strider
am 5 Jul. 2014
Bearbeitet: Star Strider
am 5 Jul. 2014
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 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!