Extracting additional parameters from a model output
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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
0 comentarios
Respuestas (1)
Cris LaPierre
el 6 de Ag. de 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)
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
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!