Solving a set of ordinary differential equations and functions dependent on each other

6 views (last 30 days)
I have these functions and ODE's from the paper “A Mathematical Model of Honey Bee Colony Dynamics to Predict the Effect of Pollen on Colony Failure” that I wanted to simulate on MATLAB, but I dont seem to know how to do that since each equation/function is really dependent on the other. I tried using ode45 to solve each on individually but this doesn't really give an accurate picture since the values aren't updating. Furthermore, I believe that I need a for loop to solve this but I have no idea how to put ode45 in a for loop or how to code the order of these equations in order to do so. Does anyone have any advice?
Bagheri, Shahin, and Mehdi Mirzaie. “A Mathematical Model of Honey Bee Colony Dynamics to Predict the Effect of Pollen on Colony Failure.” PLOS ONE, vol. 14, no. 11, 22 Nov. 2019, p. e0225632, 10.1371/journal.pone.0225632. Accessed 1 June 2020.
  3 Comments
Luis Gonzalez Barranca
Luis Gonzalez Barranca on 1 Jun 2021
very sorry about this, i had made an honest mistake writing PDE instead of ODE in the title. The code is meant to look at the effects of the survival and recruitment functions on the differential equations with defined amounts of pollen and nectar in the colony. I'm not entirely sure about the relationship between the initial condition vector and the ode solver. This has been the biggest issue so far with the following error popping up a lot:
@(T,X)[(L.*S1)-((PHI_O).*(X(1)));(PHI_O.*X(1))-(PHI_C.*(X(2)))-(MC.*(X(2)));(PHI_C.*X(2))-(X(3).*RP1)-(X(3).*RN1);(X(3).*RP1)-(MP.*X(4));(X(3).*RN1)-(MN.*X(5))] returns a vector of length 401, but the length of initial conditions vector is 5. The vector returned by @(T,X)[(L.*S1)-((PHI_O).*(X(1)));(PHI_O.*X(1))-(PHI_C.*(X(2)))-(MC.*(X(2)));(PHI_C.*X(2))-(X(3).*RP1)-(X(3).*RN1);(X(3).*RP1)-(MP.*X(4));(X(3).*RN1)-(MN.*X(5))] and the initial conditions vector must have the same number of elements.
this is the code that I have so far:
%% variables
% we are modeling how fp and fn affect the number of bees in the colony,
% therefore we can look at for example if we keep fp constant to the
% original value and fn is driven to 0 over x time where the matrix matches
% the dimensions that the solver needs
% we can look at driving fn down slowly or quickly
phi_c = 1/12; % the pupation rate of capped brood that changes to bee per day
phi_o = 1/9; % pupation rate of uncapped brood that changes to pupae per day
mc = 0.15; % look this up
mn = 0.35; % look this up
mp = 0.75; % look this up
L = 2000; % the number of eggs laid daily by the queen
v = 5000; % the number of hive bees for 50% egg survival
fn1 = 100; % nectar in grams
fp1 = 100; % pollen in grams
b = 500; % mass of nectar stored for 50% egg survival
K = 8; % the maximum protein that can saturate a nurse bee
a_minp = 0.25; % hive bee is recruited to become a pollen forager
a_maxp = 0.25; % hive bee is recruited to become a pollen forager
delta = 0.75; % effect of excess foragers on recruitment
a_minn = 0.25; % hive bee is recruited to become a nectar forager
a_maxn = 0.25; % hive bee is recruited to become a nectar forager
Fn1 = 100; % the number of nectar forager bees
Fp1 = 100; % the number of pollen forager bees
H1 = 100; % numbers of hive bees
Bo = 10000; % The number of eggs and larvae (uncapped brood)
t = linspace(1,1000);
gBo = 0.018; %daily pollen requirement per uncapped brood
c = 0.1; %food gathered per day per forager
gH = 0.0007; %daily pollen requirement per hive bee
lBo = 0.018; %daily nectar requirement per uncapped brood
lA = 0.0007; %daily nectar requirement per adult bee
fpm = linspace(0,100)'; %this is a matrix for different fp values
fnm = linspace(0,100)'; %this is a matrix for different fn values
% Function 5: The pollen recruitment function
%Rp(H, Fp, Fn, fp)
Rp1 = a_minp;
Rp2 = a_maxp*(1 - ((fpm.^2)/((fpm.^2)+(K*H1))));
Rp3 = delta * (Fp1/(Fp1+Fn1+H1));
Rp = Rp1 + Rp2 - Rp3;
% Function 6: The nectar recruitment function
%Rn(H, Fp, Fn, fn)
Rn1 = a_minn;
Rn2 = a_maxn*(1 - (fnm./(fnm+b)));
Rn3 = delta * (Fn1/(Fp1+Fn1+H1));
Rn = Rn1 + Rn2 - Rn3;
%Rn = a_minn + a_maxn*(1 - (fn1/(fn1+b))) - delta * (Fn1./(Fp1+Fn1+H1));
%survival function
%S(H, fp, fn)
S = (H1/(H1+v)).*(fnm./(fnm+b)).*((fpm.^2)/((fpm.^2)+(K*H1)));
%for S and Rp just call on the last column
S1 = S(:,100);
Rp1 = Rp(:,100);
Rn1 = Rn;
%%solving ODE's
x10 = 100;
x20 = 100;
x30 = 1000;
x40 = 100;
x50 = 100;
syms t x
ode = @(t, x) [(L.*S1) - ((phi_o).*(x(1))); ...
(phi_o.*x(1)) - (phi_c.*(x(2))) - (mc.*(x(2))); ...
(phi_c.*x(2)) - (x(3).*Rp1) - (x(3).*Rn1); ...
(x(3).*Rp1) - (mp.*x(4)); ...
(x(3).*Rn1)- (mn.*x(5))];
[t,x] = ode45(ode, [0,1000], [x10 x20 x30 x40 x50]);
Luis Gonzalez Barranca
Luis Gonzalez Barranca on 1 Jun 2021
if anyone has any suggestions for how to rework the code so that the survival and recruitment functions can update with the solutions to the ode's that would be very helpful!

Sign in to comment.

Answers (1)

Steven Lord
Steven Lord on 1 Jun 2021
Let V = [B_o; B_c; H; F_p; F_n; f_p; f_n]. I think I read those correctly from the image. Then write your function so it computes dV/dt as a function of t and V.
function dVdt = mysystem(t, V)
B_o = V(1);
B_c = V(2);
% etc
dVdt = zeros(size(V));
% Now define the elements of dVdt in terms of B_o, B_c, etc.
% You can define other intermediate variables as you need them too.
% I'd consider putting the Description from each row of the table as a
% comment before you implement that row, for easy traceability between your
% implementation and the paper's description of what you're implementing.
dVdt(1) = % dB_o/dt
end
  1 Comment
Luis Gonzalez Barranca
Luis Gonzalez Barranca on 2 Jun 2021
I'm seeing that this would help to create a place for all the differential equations but then how would i be able to solve all of them simultaneously and then update the values of my functions (second picture) with the solutions to the differential equations at the same time if i use a function?

Sign in to comment.

Categories

Find more on Programming 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!