What is the best way to change the matrix(dMdt in my code) dimension to plot it?

2 views (last 30 days)
function mu
% 6 ODEs:
% dMsh/dt = kG*Mshc*Mshn/Msh - sigmaSS*klitt*Msh/(1 + KMlitt/Msh);
% dmrt/dt = kG*Mrtc*Mrtn)/Mrt - sigmaSS*klitt*Mrt/(1+KMlitt/Mrt);
% dMshC/dt = kc*Msh/((1 + sigmaSS*Msh/KM)*(1 + sigmaPI*Mshc/(Msh*Jc)))...
% -fc*kG*Mshc*Mshn/Msh - (Mshc/Msh - Mrtc/Mrt)/(pc/Msh.^q + pc/Mrt.^q);
% dmrtC/dt = (Mshc/Msh - Mrtc/Mrt)/((pc/Msh.^q + pc/Mrt.^q) - fc*kG*Mrtc*Mrtn/Mrt;
% dMshN/dt = (Mrtn/Mrt - Mshn/Msh)/(pn/Mrt.^q + pn/Msh.^q) - fc*kG*Mshc*Mshn/Msh;
% dmrtN/dt = kn*Mrt/((1+sigmaSS*Mrt/KM) * (1+sigmaPI*Mrtn/(Mrt*Jc)))...
%- fn*kG*Mrtc*Mrtn/Mrt - (Mrtn/Mrt - Mshn/Msh)/(pn/Mrt.^q + pn/Msh.^q);
%Define Constants
kG = 200; % Growth Rate
klitt = 0.05; % Litter Rate
KMlitt = 0.5; % Litter Parameter
Jc = 0.1; % Inhibition of C assimilation
Jn = 0.01; % Inhibition of N uptake
fc = 0.5; % Fractions of C in structural
fn = 0.025; %% fN 0.025 is better % Fractions of N in structural
pc = 1; % Transport resistance coefficients
pn = 1; % Transport resistance coefficients
q = 1; % Transport resistance scaling parameter
kc = 0.1; % C assimilation parameter
kn = 0.02; % N uptake parameter
KM = 1; % Parameter giving asymtotic values of photosynthesis and N uptake
%Switches
sigmaPI_vec = [0,1]; % Product Inhibition 1(on) 0(off)
sigmaSS = 0; % Steady state growth
for i = 1:2
sigmaPI = sigmaPI_vec(i);
%Initial values defined in problem
Msh0 = 0.8;
Mrt0 = 0.9;
Mshc0 = 0;
Mrtc0 = 0.03;
Mshn0 = 0.04;
Mrtn0 = 0.001;
x0 = [Msh0 , Mrt0 , Mshc0 , Mrtc0 , Mshn0 , Mrtn0 ];
%% First solve from 0 to 400
tspan1 = [0, 400];
x01 = x0;
[t1,x1] = ode15s(@fun_call,tspan1,x01);
mass_old = x1(:,1)+x1(:,2);
%% Second solve from 400 to 500 (reduce shoot by 75%)
x02 = x1(end,:);
% Apply 75% reduction
x02(1) = 0.25*x02(1); % Msh
x02(3) = 0.25*x02(3); % MshC
x02(5) = 0.25*x02(5); % MshN
tspan2 = [400, 500];
[t2 , x2] = ode15s(@fun_call,tspan2,x02);
mass_new = x2(:,1) + x2(:,2);
%% Combine results for plotting
t_all = [t1;t2];
x_all = [x1;x2];
M = [mass_old; mass_new]; % For figure 2_Part A Specific Growth Rate
dMdt = diff(M)./diff(t_all);
%dMdt(end+1) = NaN; % Add length
mu = dMdt./M; % Specific Growth Rate of Plant
figure (2)
hold on
plot(t_all, mu, 'LineWidth', 2)
%axis([380 500 0.00 0.14])
grid on
legend('No Product Inhibition','Product Inhibition');
xlabel('Time')
ylabel('Specific Growth Rate, mu')
title('A, Specific Growth Rate')
end
function u = fun_call(~,x)
u = zeros(6,1);
Q1 = x(3)/x(1) - x(4)/x(2);
Q2 = pc/x(1)^q + pc/x(2)^q;
R1 = x(6)/x(2) - x(5)/x(1);
R2 = pn/x(2)^q + pn/x(1)^q;
S = sigmaSS*klitt;
u(1) = kG*x(3)*x(5)/x(1) - S*x(1)/(1 + KMlitt/x(1));
u(2) = kG*x(4)*x(6)/x(2) - S*x(2)/(1+KMlitt/x(2));
u(3) = kc*x(1)/((1 + sigmaSS*x(1)/KM)*(1 + sigmaPI*x(3)/(x(1)*Jc)))...
-fc*kG*x(3)*x(5)/x(1) - Q1/Q2;
u(4) = Q1/Q2 - fc*kG*x(4)*x(6)/x(2);
u(5) = R1/R2 - fn*kG*x(3)*x(5)/x(1);
u(6) = kn*x(2)/((1 + sigmaSS*x(2)/KM) * (1 + sigmaPI*x(6)/(x(2)*Jn)))...
- fn*kG*x(4)*x(6)/x(2) - R1/R2;
end
end

Answers (0)

Categories

Find more on General Applications in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!