Extracting additional parameters from a model output

1 view (last 30 days)
I have written the following code for a basic algae-phosphorous model. I want to extract the values of Pl (P limitation) to be used somewhere else. When I try to declare Pl as a global variable, it only returns a single value as opposed to an array for all the iterations. I would be very much grateful if I can have help in sorting this out.
Thanks in advance.
global mu M Ph RP SD K SR DR DD Pl
mu = 0.25; %maximum growth rate of phytoplankton (per day)
M = 0.2; % Mortality rate (per day)
Ph = 0.03; % Half saturation constant of P (umol/kg)
RP = 0.2/365;% river input (mmol per m2 per day)
SD = 500; %surface layer depth (m)
K = 3/365; %Mixing coefficient (m per day)
SR = 0.95; %shallow water P regeneration fraction
DR = 0.048; %deep water P regeneration fraction
DD = 3230; %deep layer depth (m)
V0 = [0.2; 0.1; 0.05];
tic
[t, V] = ode45(@func_p,[0 35000000],V0);
toc
plot(t,V,'-o','LineWidth',1,'MarkerSize',3)
hold on
xlim([1 35000000]);
ylabel('concentration (umol/kg)');
legend ('surface P','deep P','Algae');
xlabel('days')
grid on
function dvdt = func_p(t,V)
global mu M Ph RP SD K SR DR DD Pl;
Ps=V(1);
Pd=V(2);
A=V(3);
Pl = Ps/(Ps+Ph);
dPsdt = (RP/SD)+((K*(Pd-Ps))/SD) - (mu*Pl*A)+(M*SR*A);
dPddt = (M*DR*A*(SD/DD))-((K*(Pd-Ps))/DD) ;
dAdt = (mu*Pl-M)*A;
dvdt = [dPsdt; dPddt; dAdt];
end

Answers (1)

Cris LaPierre
Cris LaPierre on 6 Aug 2021
Global does not work because Pl only contains a single value at a time. You could try one of the solution in this answer, but I found it slow.
You have all the data you need to recalculate Pl yourself from V and Ph. Why not just do that?
Ps = V(:,1);
Pl = Ps./(Ps+Ph)
Here's an example
mu = 0.25; %maximum growth rate of phytoplankton (per day)
M = 0.2; % Mortality rate (per day)
Ph = 0.03; % Half saturation constant of P (umol/kg)
RP = 0.2/365;% river input (mmol per m2 per day)
SD = 500; %surface layer depth (m)
K = 3/365; %Mixing coefficient (m per day)
SR = 0.95; %shallow water P regeneration fraction
DR = 0.048; %deep water P regeneration fraction
DD = 3230; %deep layer depth (m)
V0 = [0.2; 0.1; 0.05];
[t, V] = ode45(@func_p,[0 35000000],V0,[],mu,M,Ph,RP,SD,K,SR,DR,DD);
Pl = V(:,1)./(V(:,1)+Ph)
Pl = 132021×1
0.8696 0.8672 0.8648 0.8622 0.8594 0.8488 0.8366 0.8235 0.8103 0.7984
plot(t,V,'-o','LineWidth',1,'MarkerSize',3)
xlim([1 35000000]);
ylabel('concentration (umol/kg)');
legend ('surface P','deep P','Algae');
xlabel('days')
grid on
function dvdt = func_p(t,V,mu,M,Ph,RP,SD,K,SR,DR,DD)
Ps=V(1);
Pd=V(2);
A=V(3);
Pl = Ps/(Ps+Ph);
dPsdt = (RP/SD)+((K*(Pd-Ps))/SD) - (mu*Pl*A)+(M*SR*A);
dPddt = (M*DR*A*(SD/DD))-((K*(Pd-Ps))/DD) ;
dAdt = (mu*Pl-M)*A;
dvdt = [dPsdt; dPddt; dAdt];
end
  1 Comment
Amavi Silva
Amavi Silva on 6 Aug 2021
Thank you! I get your point. But the thing is, using this simple model as an example, I am trying to build a more complex biogeochemical model with number of parameters influencing each other and the state variables. I want some of these paramaters out of the function and recalculating them seperately does not feel like an option due to the complexity of the model. Therefore, I thought if I can find a way to extract Pl, I can apply the same method to my complex model. Anyway, I will go through the answer you have attached here. Thank you so much again!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by