ode45 For Three Functions

Hello! I am trying to run a code that will use ode45, but I have three functions. I do not know how to write a code to combine all three. The code I have written is shown below.

I'm also hoping you can give me details as to where in a file it should be - such as if I should combine the functions into one script or what. Each of my functions is a separate function script and the main code is its' own script as well.

First function

function dHdt = hive(t,H)

% This function is used to model the population changes in the hive bees

L = 1780; % Maximum egg laying rate of queen bee
omega = 27000; % Brood mortality
alpha = 0.15; % Maximum rate at which hive bees become foragers
sigma = 0.75; % social inhibition
ro = 0.000005; % death rate due to mites

% H = population of hive bees
% F = population of forager bees
% M = population of mites

dHdt = (L*(H+F))/(omega+H+F) - H*(alpha-sigma*(F/H+F))-ro*M*H;
end

Second Function

function dFdt = forager(t,F)

% This function is used to calculate the changes in the forager bee
% population

alpha = 0.15; % Maximum rate at which hive bees become foragers
sigma = 0.75;  % Social inhibition
m =  0.12;   % Per capita deat rate of foragers
ro = 0.000005;  % Death rate due to mites

% H = population of hive bees
% F = Population of forager bees
% M = population of mites

dFdt = H*(alpha-sigma*(F/H+F))-m*F-ro*M*F;
end

Third Function

function dMdt = mite(t,M)

% This function calculates the change in the mite populations over time

r = 0.0165;    % Growth rate of mites
alphat =  0.3;  % Carrying capacity of mites in the hive
rotwo = 0.016;    % Per capita death rate of mites

% M = Population of mites
% H = Population of hive bees
% F = Population of forager bees


dMdt = r*M*(1-M/alphat*H)-rotwo*M;
end

Main Code

% What will happen to the population of bees in 12 months?

Hinitial = 14000;   % Initial Population of Hive Bees
Finitial = 7000;    % Initial population of forager bees
Minitial = 400;     % Initial population of mites
tspan = [0,12];     % Timespan over which the data will be calculated

[T, H, F, M] = ode45(@hive, @forager, @mite, tspan, Hinitial, Finitial, Minitial);

figure
plot(T,H) 
xlabel('Time (Years)')
ylabel('Number of Hive Bees')

plot(T,F) 
xlabel('Time (Years)')
ylabel('Number of Forager Bees')

plot(T,M) 
xlabel('Time (Years)')
ylabel('Number of Mites')

I need the ability to make it into a plot that shows all three populations as well.

Thank you!

Answer 1

[Alan Stevens]

Well, this is how you could structure it all in one file

% What will happen to the population of bees in 12 months? Months or years?

Hinitial = 14000;   % Initial Population of Hive Bees
Finitial = 7000;    % Initial population of forager bees
Minitial = 400;     % Initial population of mites
Binitial = [Hinitial, Finitial, Minitial];
tspan = [0 12];     % Timespan over which the data will be calculated

[T, B] = ode45(@bees,tspan,Binitial);
H = B(:,1);
F = B(:,2);
M = B(:,3);

figure
plot(T,H) 
xlabel('Time (Years)')  %% Years or months? 
ylabel('Number of Hive Bees')
figure
plot(T,F) 
xlabel('Time (Years)')
ylabel('Number of Forager Bees')
figure
plot(T,M) 
xlabel('Time (Years)')
ylabel('Number of Mites')

function dBdt = bees(~,B)

% This function is used to model the population changes 

L = 1780; % Maximum egg laying rate of queen bee
omega = 27000; % Brood mortality
alpha = 0.15; % Maximum rate at which hive bees become foragers
sigma = 0.75; % social inhibition
ro = 0.000005; % death rate due to mites
r = 0.0165;    % Growth rate of mites
alphat =  0.3;  % Carrying capacity of mites in the hive
rotwo = 0.016;    % Per capita death rate of mites
m =  0.12;   % Per capita deat rate of foragers

% H = population of hive bees
% F = population of forager bees
% M = population of mites
H = B(1);
F = B(2);
M = B(3);

dHdt = (L*(H+F))/(omega+H+F) - H*(alpha-sigma*(F/H+F))-ro*M*H;
dFdt = H*(alpha-sigma*(F/H+F))-m*F-ro*M*F;
dMdt = r*M*(1-M/alphat*H)-rotwo*M;

dBdt = [dHdt; dFdt; dMdt];
end

However, the results don't look sensible! You should check your equations carefully, as well as making sure that the units of your constants (rates) are all consistent.

[Anneliese Fensch]

Thank you for helping me out! It does seem that there is something up with my equations. I'm trying to model something from a paper I read for class. I am probably typing them wrong.

I really was mainly concerned about how to get ode45 to run multiple functions so I can figure that out in the future.

Thank you so much!

[Anneliese Fensch]

If you don't mind, I'm going to show you the equations that I am trying to type. Could you help me with that too?

[Alan Stevens]

Yes, there were some errors in your implementation of the equations. Try

% What will happen to the population of bees in 12 months? Months or years?

Hinitial = 14000;   % Initial Population of Hive Bees
Finitial = 7000;    % Initial population of forager bees
Minitial = 400;     % Initial population of mites
Binitial = [Hinitial, Finitial, Minitial];
tspan = [0 12];     % Timespan over which the data will be calculated

[T, B] = ode45(@bees,tspan,Binitial);
H = B(:,1);
F = B(:,2);
M = B(:,3);

figure
plot(T,H),grid
xlabel('Time (Years)')  %% Years or months? 
ylabel('Number of Hive Bees')
figure
plot(T,F),grid
xlabel('Time (Years)')
ylabel('Number of Forager Bees')
figure
plot(T,M),grid 
xlabel('Time (Years)')
ylabel('Number of Mites')

function dBdt = bees(~,B)

% This function is used to model the population changes 

L = 1780; % Maximum egg laying rate of queen bee
omega = 27000; % Brood mortality
alpha = 0.15; % Maximum rate at which hive bees become foragers
sigma = 0.75; % social inhibition
ro = 0.000005; % death rate due to mites
r = 0.0165;    % Growth rate of mites
alphat =  0.3;  % Carrying capacity of mites in the hive
rotwo = 0.016;    % Per capita death rate of mites
m =  0.12;   % Per capita deat rate of foragers

% H = population of hive bees
% F = population of forager bees
% M = population of mites
H = B(1);
F = B(2);
M = B(3);

dHdt = (L*(H+F))/(omega+H+F) - H*(alpha-sigma*(F/(H+F)))-ro*M*H;
dFdt = H*(alpha-sigma*(F/(H+F)))-m*F-ro*M*F;
dMdt = r*M*(1-M/(alphat*H))-rotwo*M;

dBdt = [dHdt; dFdt; dMdt];
end

function dFdt = forager(~,F)

% This function is used to calculate the changes in the forager bee
% population

alpha = 0.15; % Maximum rate at which hive bees become foragers
sigma = 0.75;  % Social inhibition
m =  0.12;   % Per capita deat rate of foragers
ro = 0.000005;  % Death rate due to mites

% H = population of hive bees
% F = Population of forager bees
% M = population of mites

dFdt = H*(alpha-sigma*(F/H+F))-m*F-ro*M*F;
end

function dMdt = mite(~,M)

% This function calculates the change in the mite populations over time

r = 0.0165;    % Growth rate of mites
alphat =  0.3;  % Carrying capacity of mites in the hive
rotwo = 0.016;    % Per capita death rate of mites

% M = Population of mites
% H = Population of hive bees
% F = Population of forager bees


dMdt = r*M*(1-M/alphat*H)-rotwo*M;
end

[Anneliese Fensch]

Thank you so much! That got me the exact graphs that I needed.

Answer 2

[James Tursa]

You made a good start, but you need to solve these differential equations simultaneously. That will mean you need to work with a state vector. Using the nomenclature of the ode45 doc, call this state vector y. Since you have three variables of interest (H, F, M), that will mean a 3-element state vector. So define

y(1) = H

y(2) = F

y(3) = M

Then define your derivative function in terms of y. You could combine all three of your current functions into one function, or call all three functions individually and combine the results. To keep changes to your current code at a minimum since you have already written the three individual derivative functions, I will show you the latter approach. So,

Step 1:

Change your function signatures to pass in H, F, and M from y to use in your derivative code since all three functions depend on all three of these variables. E.g.,

function dHdt = hive(t,H,F,M)
etc.
function dFdt = forager(t,H,F,M)
etc.
function dMdt = mite(t,H,F,M)
etc.

Step 2:

In your main code, create a function handle with a (t,y) signature that passes in the H, F, M from the 3-element y vector and returns a 3-element derivative vector. E.g., since we have defined y(1)=H, y(2)=F, and y(3)=M, this would result in:

f = @(t,y) [hive(t,y(1),y(2),y(3));forager(t,y(1),y(2),y(3));mite(t,y(1),y(2),y(3))]

Step 3:

In your main code, call ode45 with the f function handle and with a vector initial condition. E.g.,

[T, Y] = ode45(f, tspan, [Hinitial;Finitial;Minitial]);

Then pick off the appropriate columns of Y for the H, F, and M solutions.