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

1 view (last 30 days)
AT on 12 May 2022
Commented: AT on 12 May 2022
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
AT on 12 May 2022
I have attached the correct plot part. Mine is not like this.

### Categories

Find more on General Applications in Help Center and File Exchange

R2020b

### Community Treasure Hunt

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

Start Hunting!