How to store the results of a system of odes solved in a loop?

2 次查看（过去 30 天）
Charlotte Coates
Charlotte Coates 2016-11-20
评论： Star Strider 2016-11-23
I'm trying to use the results of the ODE Solver to calculate a variable's new value and store all of the results so that the variables can be plotted against time. The variables I am calculating(k, Q) are functions of the solutions to the system of differential equations and I want to save them as vectors. I'm not sure how to use the solution of each iteration of loop.
%---- Initial conditions
y0 =[10; 10]; % [Initial infected tillers; Initial uninfected tillers]
k(1) = 1 - ((10+10)/4000); % Scalar of maxiumum rates
Q(1) = y0(1)/(y0(1)+y0(2)); % Proportion of infected tillers
for i = 2 : 8
a_b(i) = a_b(i-1) + 0.1; % Varying birth advantages
deq1=@(t,y)[(1 - lambda)*k.*a_b(i).*b.*y(1) - y(1).*(a_d.*d); k.*b.*y(2)+lambda*k.*a_b(i).*y(1) - d.*y(2)];
[t,sol] = ode15s(deq1, t_int, y0);
k(i) = (1 - (sol(:,1)+ sol(:,2)))/4000; % Scalar of maxiumum growth rates
arraysize = size(t);
Q(i) = (sol(:,1))/(sol(:,1)+ (sol(:,2))); % Find the quasi equilibrium
end

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

回答（1 个）

Star Strider
Star Strider 2016-11-20
If you want to save the ‘t’ and ‘sol’ in each loop iteration, use a cell array:
[t{i},sol{i}] = ode15s(deq1, t_int, y0);
  9 个评论
显示 7更早的评论
Charlotte Coates
Charlotte Coates 2016-11-23
I actually need to calculate k and Q for every value of t, not for every value of i. So they are vectors of t = size(t) magnitude. How can use the output of ode15s for each value of t to determine the k and Q values so that the following integration uses the new k value for time (t+0.25) I calculate with the populations at time t.
Star Strider
Star Strider 2016-11-23
My pleasure!
One possible problem I notice is that you subscript ‘k’, but you also use ‘k’ in unsubscripted form in ‘deq1’:
[(1 - lambda)*k.*a_b(i) ...
Note that ‘k(i)’ is a scalar and ‘k’ is a vector. The vector may not throw an error if everything else works, but may not give you the result you want. Which one (of ‘k’ or ‘k(i-1)’) do you want to use in ‘deq1’?
I’m not certain what you’re modeling, so I can’t claim any expertise on your code. One other possibility is that there’s a subtle problem in your code, perhaps with the value of a constant. If other people have modeled this system, see if they encountered negative values that shouldn’t be negative, and how they dealt with that. Run your simulations with their data to see if you get the same results they did.
I notice that ‘deq1’ returns a (2x1) vector. If the first row is a derivative of the second row (I have no idea), then you can expect the derivative to go negative if the function itself decreases.
It might be possible to ‘trap’ a negative value and reset it, or use the absolute value or zero instead. The very real problem with this approach is that it introduces discontinuities, and the numeric solvers have serious problems with discontinuities.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by