Can someone check my code and rectify me, getting constant output but needs changing.

1 view (last 30 days)
Can someone check my code and rectify my error.
Here I have given a main code, two functions. On running it gives me constant output but I need that value of 'a' in main function should keep changing with each loop and final value I need like Glc get stored, at end of loop, to plot a variable graph. I have tried to simplify code and added possible change in comments format, but there's still problem. So, I needed help in it. In advance Thanks a lot for helping me out.
clear
close all
clc
M = {'Gal','Urd','Glc','Asn','Amm','Glu','Lac','Gln','Asp'};
v = 100;
cellden = 200000;
mb = 2.76709;
umax = 0.0650;
gal = 0;
urd = 0;
glc = 36;
asn = 6.09;
amm = 0.1;
glu = 2.125333333;
lac = 1.391166667;
gln = 0.01;
asp = 1.05;
for k = 1:216
a = flim(glc,asn);
b = 2;
u = umax*a*b;
udeathf = 3;
qAmm = 12;
qUrd = 19;
qGln = 15;
qGlc = -10;
qGal = 11;
qLac = 12;
qAsn = -13;
qAsp = -12;
qGlu = -9;;
tspan = [0 216];
xo = [100;20000000;0;0;3600;609;10;212.533;139.11667;1;105;276.709];
[t,X]=ode23s(@(t,X) new_try_func(t,X,v,cellden,mb,u,udeathf,gal,urd,glc,asn,amm,glu,lac,gln,asp,qAmm,qUrd,qGln,qGlc,qGal,qLac,qAsn,qAsp,qGlu),tspan,xo);
cd(k) = X(1,2)/X(1,1);
time(k) = k;
v = X(1,1);
cellden = X(1,2)/X(1,1);
gal = X(1,3)/X(1,1);
urd = X(1,4)/X(1,1);
glc = X(1,5)/X(1,1);
asn = X(1,6)/X(1,1);
amm = X(1,7)/X(1,1);
glu = X(1,8)/X(1,1);
lac = X(1,9)/X(1,1);
gln = X(1,10)/X(1,1);
asp = X(1,11)/X(1,1);
mb = X(1,12)/X(1,1);
end
%cell_density(:,1) = X(:,2)./X(:,1);
%Gal(:,1) = X(:,3)./X(:,1);
%Urd(:,1) = X(:,4)./X(:,1);
%Glc(:,1) = X(:,5)./X(:,1);
%Asn(:,1) = X(:,6)./X(:,1);
%Amm(:,1) = X(:,7)./X(:,1);
%Glu(:,1) = X(:,8)./X(:,1);
%Lac(:,1) = X(:,9)./X(:,1);
%Gln(:,1) = X(:,10)./X(:,1);
%Asp(:,1) = X(:,11)./X(:,1);
%mAb(:,1) = X(:,12)./X(:,1);
plot(time,cd); % or plot(cell_density);
xlim([0 216]);
Next are function code in simplified way(not original)
function dydt = new_try_func(~,X,v,cellden,mb,u,udeathf,gal,urd,glc,asn,amm,glu,lac,gln,asp,qAmm,qUrd,qGln,qGlc,qGal,qLac,qAsn,qAsp,qGlu)
mmAb = 0.410;
YmAb_X = 3.39;
dydt = zeros(12,1);
Fin = 0.1;
Fout = 0;
Gal = 0;
Urd = 0;
Glc = 144.37;
Asn = 26.99;
Amm = 0.06;
Glu = 12.19;
Lac = 16.64;
Gln = 0;
Asp = 51.95;
dydt(1,1) = Fin - Fout; % volume
dydt(2,1) = ((u - udeathf)*v*cellden) - (Fout*cellden); %cell density
dydt(3,1) = (Fin*Gal) - (Fout*gal) + (qGal*v*cellden); % Gal
dydt(4,1) = (Fin*Urd) - (Fout*urd) + (qUrd*v*cellden); % Urd
dydt(5,1) = (Fin*Glc) - (Fout*glc) + (qGlc*v*cellden); % Glc
dydt(6,1) = (Fin*Asn) - (Fout*asn) + (qAsn*v*cellden); % Asn
dydt(7,1) = (Fin*Amm) - (Fout*amm) + (qAmm*v*cellden); % Amm
dydt(8,1) = (Fin*Glu) - (Fout*glu) + (qGlu*v*cellden); % Glu
dydt(9,1) = (Fin*Lac) - (Fout*lac) + (qLac*v*cellden); % Lac
dydt(10,1) = (Fin*Gln) - (Fout*gln) + (qGln*v*cellden); % Gln
dydt(11,1) = (Fin*Asp) - (Fout*asp) + (qAsp*v*cellden); % Asp
dydt(12,1) = (-1*Fout*mb) + ((YmAb_X*u + mmAb)*v*cellden); %mAb
%v = X(1);
%cellden = X(2)/X(1);
%gal = X(3)/X(1);
%urd = X(4)/X(1);
%glc = X(5)/X(1);
%asn = X(6)/X(1);
%amm = X(7)/X(1);
%glu = X(8)/X(1);
%lac = X(9)/X(1);
%gln = X(10)/X(1);
%asp = X(11)/X(1);
%mb = X(12)/X(1);
end
Function for 'a'
function a = flim(glc,asn)
KGlc = 14;
KAsn = 2;
a = ((glc*asn)/((glc + KGlc)*(asn + KAsn)));
end

Accepted Answer

Walter Roberson
Walter Roberson on 12 Jun 2021
Edited: Walter Roberson on 12 Jun 2021
[t,X]=ode23s(@(t,X) new_try_func(t,X,v,cellden,mb,u,udeathf,gal,urd,glc,asn,amm,glu,lac,gln,asp,qAmm,qUrd,qGln,qGlc,qGal,qLac,qAsn,qAsp,qGlu),tspan,xo);
There is one row of X for each t value. The rows correspond to the X vector at different time steps.
glc = X(1,5)/X(1,1);
The different rows correspond to the X vector at different time steps. The first row, X(1,:), corresponds to the X values at the first time step. Which is to say, it is always going to be the same as the boundary conditions you pass in.
You are therefore going to be setting the new value of glc the same as what you passed in for the boundary conditions. You are passing in 100 for x0(1) and 3600 for x0(6) so you will get out 3600/100 = 36 which is the same as your initial glc. Nothing can change.
You should be considering extracting from (for example) X(end,5)/X(end,1)
(My internal tracking number for this is T0098912)
  2 Comments
Agrim Singh
Agrim Singh on 12 Jun 2021
Edited: Agrim Singh on 12 Jun 2021
The values are updating now, BUT IT stops updating after t = 24; and becomes constant. Any clues to why???
Walter Roberson
Walter Roberson on 14 Jun 2021
glc and asn get up to roughly 1e154, and when that happens their product is infinite and the ratio of two infinite numbers in flim leads to NaN.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!