Extracting additional parameters from a model output

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

Answer 1

[Cris LaPierre]

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

[Amavi Silva]

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!