Dear all.
I am trying to use the ndsolve function to solve a system of ODEs, but I can't since I always obtain "singularity/stiffness system" errors.
I would be very glad if someone can advice me or give me some assistance about this issue.
For more detail, I attach the Mathematica code and the obtained errors are summarized next:
NDSolve::ndsz : At t == 0.7443051546547935, step size is effectively zero; singularity or stiff system suspected.
InterpolatingFunction::dmval : Input value {0.8} lies outside the range of data in the interpolating function. Extrapolation will be used.
InterpolatingFunction::dmval : Input value {0.9} lies outside the range of data in the interpolating function. Extrapolation will be used.
InterpolatingFunction::dmval : Input value {1.} lies outside the range of data in the interpolating function. Extrapolation will be used.
General::stop : Further output of InterpolatingFunction::dmval will be suppressed during this calculation.
Thanks a lot.
Best regards.
Manuel.
Mathematica code:
raftot = 3/1000;(*.003;*)
prmt5 = 3/1000;(*.003;*) (*0.0*)
e2tot = 3/10000;(*.0003;*)
mektot = 12/10;(*1.2;*)
mekpasetot = 3/10000;(*.0003;*)
erktot = 12/10;(*1.2;*)
erkpasetot = 12/100;(*.12;*)
a1 = 958.64;(*1000*)
a2 = 958.64;(*1000*)
a3 = 958.64;(*1000*)
a4 = 958.64;(*1000*)
a5 = 958.64;(*1000*)
a6 = 958.64;(*1000*)
a7 = 958.64;(*1000*)
a8 = 958.64;(*1000*)
a9 = 958.64;(*1000*)
a10 = 958.64;(*1000*)
aprmt5 = 482.94;
d1 = 150;
d2 = 150;
d3 = 150;
d4 = 150;
d5 = 150;
d6 = 150;
d7 = 150;
d8 = 150;
d9 = 150;
d10 = 150;
k1 = 150;
k2 = 150;
k3 = 150;
k4 = 150;
k5 = 150;
k6 = 150;
k7 = 150;
k8 = 150;
k9 = 150;
k10 = 150;
i = -1.5;
(*For[i=-6,i\[LessEqual]-6,i=i+diff,*)
e1tot = 10^(i);
alphaprmt5 = e1tot/10^(-1); (*alphaprmt5=((e1tot-mose1[t])/10^(-1)*)
s = NDSolve[{
raf'[t] == -a1*raf[t]*e1[t] + d1*rafe1[t] + k2*rafpe2[t],
rafe1'[t] == a1*raf[t]*e1[t] - (d1 + k1) rafe1[t],
rafp'[t] == -a2*rafp[t]*e2[t] + d2*rafpe2[t] +
k1*rafe1[t] + (k3 + d3)*mekrafp[t] -
a3*rafp[t]*mek[t] + (k5 + d5)*mekprafp[t] -
a5*mekp[t]*rafp[t],(*-aprmt5*alphaprmt5*prmt5*mosstar[t],*)
rafpe2'[t] == a2*rafp[t]*e2[t] - (d2 + k2)*rafpe2[t],
mek'[t] == -a3*mek[t]*rafp[t] + d3*mekrafp[t] +
k4*mekpmekpase[t],
mekrafp'[t] == a3*mek[t]*rafp[t] - (d3 + k3)*mekrafp[t],
mekp'[t] == -a4*mekp[t]*mekpase[t] + d4*mekpmekpase[t] +
k3*mekrafp[t] + k6*mekppmekpase[t] + d5*mekprafp[t] -
a5*mekp[t]*rafp[t],
mekpmekpase'[t] ==
a4*mekp[t]*mekpase[t] - (d4 + k4)*mekpmekpase[t],
mekprafp'[t] == a5*mekp[t]*rafp[t] - (d5 + k5)*mekprafp[t],
mekpp'[t] ==
k5*mekprafp[t] - a6*mekpp[t]*mekpase[t] + d6*mekppmekpase[t] -
a7*mekpp[t]*erk[t] + (d7 + k7)*erkmekpp[t] + (d9 + k9)*
erkpmekpp[t] - a9*erkp[t]*mekpp[t],
mekppmekpase'[t] ==
a6*mekpp[t]*mekpase[t] - (d6 + k6)*mekppmekpase[t],
erk'[t] == -a7*erk[t]*mekpp[t] + d7*erkmekpp[t] +
k8*erkperkpase[t],
erkmekpp'[t] == a7*erk[t]*mekpp[t] - (d7 + k7)*erkmekpp[t],
erkp'[t] ==
k7*erkmekpp[t] - a8*erkp[t]*erkpase[t] + d8*erkperkpase[t] -
a9*erkp[t]*mekpp[t] + d9*erkpmekpp[t] + k10*erkpperkpase[t],
erkpmekpp'[t] == a9*erkp[t]*mekpp[t] - (d9 + k9)*erkpmekpp[t],
erkpp'[t] == -a10*erkpp[t]*erkpase[t] + d10*erkpperkpase[t] +
k9*erkpmekpp[t],
erkperkpase'[t] ==
a8*erkp[t]*erkpase[t] - (d8 + k8)*erkperkpase[t],
erkpperkpase'[t] ==
a10*erkpp[t]*erkpase[t] - (d10 + k10)*erkpperkpase[t],
e1'[t] == -a1*raf[t]*e1[t] + (d1 + k1)*rafe1[t],
e2'[t] == -a2*rafp[t]*e2[t] + (d2 + k2)*rafpe2[t],
mekpase'[
t] == -a4*mekp[t]*mekpase[t] + (d4 + k4)*mekpmekpase'[t] -
a6*mekpp[t]*mekpase[t] + (d6 + k6)*mekppmekpase[t],
erkpase'[
t] == -a8*erkp[t]*erkpase[t] + (d8 + k8)*erkperkpase'[t] -
a10*erkpp[t]*erkpase[t] + (d10 + k10)*erkpperkpase[t],
raf[0] == raftot, rafe1[0] == 0, rafp[0] == 0, rafpe2[0] == 0,
mek[0] == mektot, mekrafp[0] == 0, mekp[0] == 0,
mekpmekpase[0] == 0, mekprafp[0] == 0, mekpp[0] == 0,
mekppmekpase[0] == 0, mekpase[0] == mekpasetot, erk[0] == erktot,
erkmekpp[0] == 0, erkp[0] == 0, erkpmekpp[0] == 0, erkpp[0] == 0,
erkperkpase[0] == 0, erkpperkpase[0] == 0,
erkpase[0] == erkpasetot, e1[0] == e1tot, e2[0] == e2tot},
{raf, rafe1, rafp, rafpe2, mek, mekrafp, mekp, mekprafp,
mekpmekpase, mekpp, mekppmekpase, mekpase, erk, erkmekpp, erkp,
erkpmekpp, erkpp, erkperkpase, erkpperkpase, erkpase, e1, e2},
{t, 0, 500}, MaxSteps -> 10000,
Method -> {StiffnessSwitching}];
Plot[raf[t] /. s, {t, 0, 0.8}, PlotRange -> All]
Plot[rafp[t] /. s, {t, 0, 0.8}, PlotRange -> All]
Plot[mekpp[t] /. s, {t, 0, 0.8}, PlotRange -> All]
Plot[erkpp[t] /. s, {t, 0, 0.8}, PlotRange -> All]