Extracting additional parameters from a model output
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Cris LaPierre
2021-8-6
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!