Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I´m modelling a biological reaction thorugh a set of differential equations.

Modelo = NDSolve[{{dX[t]/dt == muapp*X[t] - 0.005*X[t], X[0] == Xo2},
{dSg[t]/dt == -Qs + 0.005*(100 - Sg[t]), Sg[0] == Sgo2},
{dAce[t]/dt == Qace - 0.005*Ace[t], Ace[0] == Aceo2},
{dBut[t]/dt == Qbut - 0.005*But[t], But[0] == Buto2},
{dAct[t]/dt == Qact - 0.005*Act[t], Act[0] == Acto2},
{dEt[t]/dt == Qet - 0.005*Et[t], Et[0] == Eto2},
{dSb[t]/dt == Qsb + 0.005*(2 - Sb[t]), Sb[0] == Sbo2}}, {X, Sg, Ace, But, Act, Et, Sb}, {t, 0, 400}, Method ->DOPRI];

Using NDSolve the system is solve with no problems. however, since this is a biological reaction, none of the solutions must be >0. I have already adjusted the parameters so that this condition is met and it works fine in matlab and modelmaker, but when solving it with mathematica X goes >0 between 200 < t< 400. In matlab I use ode15s, because is stiff at some points, but it works just as well with ode45 and ode23s, so i sat NDsolve parameters to those that mathematica tutorials say to be the same as ode45 in matlab:

DOPRIamat = {{1/5}, {3/40,9/40}, {44/45,-56/15,32/9}, {19372/6561,-25360/2187,64448/6561,-212/729}, {9017/3168,-355/33,46732/5247,49/176,-5103/18656}, {35/384,0,500/1113,125/192,-2187/6784,11/84}};
DOPRIbvec = {35/384, 0, 500/1113, 125/192,-2187/6784, 11/84, 0};
DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1};
DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200, 22/525, -1/40};
DOPRICoefficients[5, p_] := N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];
DOPRI = {"ExplicitRungeKutta", "DifferenceOrder" -> 5, "Coefficients" -> DOPRICoefficients, "StiffnessTest" -> False};

But still i get x<0. Is it possible to modifie NDSolve to work as ode15s in matlab? or can it at least be modified to keep it from giving negative values? thanks for your help, i leave down the script i use in mathematica.

Mathematica

Xo2 = 1.0605;
Sgo2 = 43.9316;
Sbo2 = 3.5023;
Aceo2 = 0.3437;
Buto2 = 2.5299;
Acto2 = 1.0589;
Eto2 = 0.2334;
ksg = 154.060527452649;
kibut = 7.00180294979766;
ksgb = 13.6974233244128;
ksb = 1.10465168560025;
ksgact = 36.1766330476263;
ksget = 2.48715131190954;
ksgsb = 0.500075731973805;
ksba = 0.398637629396357;
kaa = 13.8244477511478;
kd = 0.00478170792619247;
m = 1.04817713477516;
maxx = 0.213696652849562;
maxbut = 0.150892601110339;
maxact = 0.778875361557128;
maxet = 0.0329313766102578;
maxsb = 0.131139790008708;
maxace = 0.142882096617404;
Ybutx = 1.22938023298761;
Yactx = 0.129135371286816;
Yetx = 0.0759697564493146;
Ysbx = 1.01875350101917;
Yacex = 2.10990339926643;
Yxsg = 0.284627874526466;
Ybutsg = 0.819985406580864;
Yactsg = 0.851780607333819;
Yetsg = 0.598346992105451;
Ysbsg = 0.995215479964112;
Yaceact = 0.973079461528239;
Ybutsb = 0.839985001333182;
Yacesb = 0.525389668809545;
mux = maxx*(Sg[t]/(ksg + Sg[t]))*(1 - (But[t]/kibut)^m);
muapp = mux - kd;
mubut = maxbut*(Sg[t]/(ksgb + Sg[t]))*(Sb[t]/(ksb + Sb[t]));
muact = maxact*(Sg[t]/(ksgact + Sg[t]));
muet = maxet*(Sg[t]/(ksget + Sg[t]));
musb = maxsb*(Sg[t]/(ksgsb + Sg[t]));
muace = maxace*(Sb[t]/(ksba + Sb[t]))*(Act[t]/(kaa + Act[t]));
Qs = X[t]*((mux/Yxsg) + ((mubut*Ybutx)/Ybutsg) + ((muact*Yactx)/Yactsg) + ((muet*Yetx)/Yetsg) + ((musb*Ysbx)/Ysbsg));
Qace = muace*X[t]*Yacex;
Qbut = mubut*X[t]*Ybutx;
Qet = muet*X[t]*Yetx;
Qact = X[t]*(muact*Yactx - ((muace*Yacex)/Yaceact));
Qsb = X[t]*(musb*Ysbx - ((mubut*Ybutx)/Ybutsb) - ((muace*Yacex)/Yacesb));
share|improve this question
2  
Why do you write dX[t]/dt? Do you mean X'[t]? Same for the other diff. eqs. (Your code does not work as is.) –  Michael E2 Jul 31 at 2:04

1 Answer 1

I don't see any negative X[t] when I NDSolve this. When I Plot[X[t]/.Modelo,{t,0,400}], it LOOKS negative, but that's because Mathematica set the y-axis range from 0.15 to 1.4. If you try a PlotRange->{0,All} in there, it looks better. Maybe this is it?

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.