# Problems with NDSolve and "singularity/stiffness system" errors.

Posted 9 years ago
6501 Views
|
2 Replies
|
0 Total Likes
|
 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]
2 Replies
Sort By:
Posted 9 years ago
 Thanks Michael.Yes, I agree with you. In fact, I think that the problem is due to a deep stiff in some variables in a very short range of time. But I don't know how manage it.
Posted 9 years ago
 Could it be a singularity? Increasing WorkingPrecision gives the same result: mekpase[t] seems to go to negative infinity: (-6.9 x 10^10) at the stopping point for WorkingPrecision -> \$MachinePrecision, -3.8 x 10^24 for WorkingPrecision -> 30. I don't think you can get around a singularity. Also mekp, mekpmekpase seem to get very large. Maybe others, too.