# NDSolve gives wrong results for “stiff system”

I have a physical problem in which I want to solve for c[e] as a function of e, where c is expected to be in the range from 0 to 1. I have a ODE to solve, which is c'[e] == f[c[e], e], and f[c[e], e] is given by

(2.*10^-6 (-0.706288 - 3.24379*10^-9 e +
c[e] (0.767155 +
6.48758*10^-9 e + (-0.27 - 2.2833*10^-9 e) c[
e])))/((1.42066 - 1. c[e])^2 ((
e (-6.4542*10^-7 +
e (2.72905*10^-15 - 1.65436*10^-30 c[e])))/(-1.42066 +
1. c[e])^3 + (
1.11743*10^-30 e^2 (6.18655*10^8 + 1.42066 e +
c[e] (-6.71971*10^8 -
2.84132 e + (2.365*10^8 + 1. e) c[e]))^2)/(1.42066 -
1. c[e])^4 + (1. + (1/((1.42066 - 1. c[e])^2))
e (0.000010341 + 2.37466*10^-14 e +
c[e] (-0.0000112321 -
4.74932*10^-14 e + (3.95315*10^-6 + 1.67152*10^-14 e) c[
e])))^2))


I use NDSolve,

g = NDSolve[{c'[e] == f[c[e], e], c[0] == 0}, c, {e, -150000, 0}]


but get a message

*NDSolve::ndsz: "At e == -120259., step size is effectively zero; singularity or stiff system suspected."

The result is an interpolating function. For e > -120258, the result is good.

For example

Plot[Evaluate[c[e] /. g], {e, -120258, 0}]


gives

But if e = -120259, c instantly becomes order 10^15, which I know must be wrong.

I also plot f[c, e] (as a surface, for all possible c and e) and f[c[e], e] (as a single curve) in 3D plot.

Show[
ParametricPlot3D[{c, e, f[c, e]}, {c, 0, 1}, {e, -150000, 0},
PlotRange -> {{0, 1}, {-150000, 0}, {-0.001, 0.001}}],
ParametricPlot3D[{c[e] /. g[[1]], e, f[c[e] /. g[[1]], e]}, {e, -120259, 0},
PlotStyle -> Red,
PlotRange -> {{0, 1}, {-150000, 0}, {-0.002, 0.002}}],
BoxRatios -> {1, 1, 1}]


gives

I see the surface is mostly smooth, but at a strip region the derivative values are quite unstable as in a "stiff system". When the solution curve enters that region it goes wrong. In this case, how can I get the correct result?

Update: as suggested by @Nasser, source codes are provided below

s11t = 7.4*10^(-12);
s33t = 13.1*10^(-12);
s12t = -1.4*10^(-12);
s13t = -4.4*10^(-12);
s44t = 16.4*10^(-12);
s66t = 7.6*10^(-12);
k11t = 4400*8.85*10^(-12);
k33t = 129*8.85*10^(-12);
d15t = 564*10^(-12);
d31t = -33.4*10^(-12);
d33t = 90*10^(-12);
M0 = {
{s11t, s12t, s13t, 0, 0, 0, 0, 0, d31t},
{s12t, s11t, s13t, 0, 0, 0, 0, 0, d31t},
{s13t, s13t, s33t, 0, 0, 0, 0, 0, d33t},
{0, 0, 0, s44t, 0, 0, 0, d15t, 0},
{0, 0, 0, 0, s44t, 0, d15t, 0, 0},
{0, 0, 0, 0, 0, s66t, 0, 0, 0},
{0, 0, 0, 0, d15t, 0, k11t, 0, 0},
{0, 0, 0, d15t, 0, 0, 0, k11t, 0},
{d31t, d31t, d33t, 0, 0, 0, 0, 0, k33t}
};
M1 = {
{s11t, s12t, s13t, 0, 0, 0, 0, 0, -d31t},
{s12t, s11t, s13t, 0, 0, 0, 0, 0, -d31t},
{s13t, s13t, s33t, 0, 0, 0, 0, 0, -d33t},
{0, 0, 0, s44t, 0, 0, 0, -d15t, 0},
{0, 0, 0, 0, s44t, 0, -d15t, 0, 0},
{0, 0, 0, 0, 0, s66t, 0, 0, 0},
{0, 0, 0, 0, -d15t, 0, k11t, 0, 0},
{0, 0, 0, -d15t, 0, 0, 0, k11t, 0},
{-d31t, -d31t, -d33t, 0, 0, 0, 0, 0, k33t}
};
I0 = IdentityMatrix[9];
L0 = Inverse[M0];
L1 = Inverse[M1];
s1111 = 1;
s1122 = L0[[1, 2]]/L0[[1, 1]];
s1133 = L0[[1, 3]]/L0[[1, 1]];
s1212 = 1/2;
s1313 = 1/2;
s1341 = -1/2*L0[[7, 5]]/(L0[[4, 4]]*M0[[7, 7]] + (L0[[7, 5]])^2);
s4242 = 1;
s4343 = 1;
S = {
{s1111, s1122, s1133, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 2 s1313, 0, 2 s1341, 0, 0},
{0, 0, 0, 0, 0, 2 s1212, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, s4242, 0},
{0, 0, 0, 0, 0, 0, 0, 0, s4343}
};
B1 = Inverse[I0 + L0.(I0 - S).(M1 - M0)];
Q[c_] := Inverse[c*B1 + (1 - c)*I0];
A[c_] := (M1 - M0).B1.Q[c];
Bp1[c_] := I0 - (1 - c)*(M1 - M0).B1.Q[c].L0.(I0 - S)
Bp2[c_] := B1.Q[c];
Bp3[c_] := Transpose[Bp2[c]];
B[c_] := Bp1[c] + Bp3[c];
C2[c_] := -(1 - c)*B1.Q[c].L0.(I0 - S);
X[e_] := {0, 0, 0, 0, 0, 0, 0, 0, e};
Yds = {0, 0, 0, 0, 0, 0, 0, 0, -0.54};
Bf2[c_] := B1.Q[c].(B1 - I0).Q[c];
Bf3[c_] := Transpose[Bf2[c]];
Af[c_] := A[c] - c*(M1 - M0).B1.Q[c].(B1 - I0).Q[c];
Bf[c_] := B[c] + c*(-Bf3[c] + (M1 - M0).B1.Q[c].(I0 + (1 - c)*(B1 - I0).Q[c]).L0.(I0 - S));
Cf[c_] := C2[c] + c*(B1.Q[c].(I0 + (1 - c)*(B1 - I0).Q[c])).L0.(I0 - S);
fdriv[c_, e_] := 1/2 (X[e].Af[c].X[e] + X[e].Bf[c].Yds + Yds.Cf[c].Yds);
n = 2.0*10^(-6);
r1 = 6.83*10^4;
r2 = 1.08*10^6;
delpc[c_, e_] := D[fdriv[x, y], x] /. {x -> c, y -> e};
delpE[c_, e_] := D[fdriv[x, y], y] /. {x -> c, y -> e};
f[c_, e_] := FullSimplify[n*delpE[c, e]/((1 - fdriv[c, e]/r1)^2 + (fdriv[c, e]/r2)^2-n*delpc[c, e])];
g = NDSolve[{c'[e] == f[c[e], e], c[0] == 0}, c, {e, -150000, 0}]

-
can you post complete self contained code that one can copy and run? You say  and f[c[e], e] is given by then give an expression. Why not post the actual code itself? i.e. the definition as you have it in your notebook. –  Nasser Jun 12 '14 at 11:24
Certainly, I have update the thread. Hope this will help, thanks. –  riorita Jun 12 '14 at 19:17
Ijust tried it with g = NDSolve[{c'[e] == f[c[e], e], c[0] == 0}, c, {e, -150000, 0}, Method -> {StiffnessSwitching, Method -> {ExplicitRungeKutta, Automatic}}, AccuracyGoal -> 5, PrecisionGoal -> 4] and get no error. You might try this and see if the plot looks ok now –  Nasser Jun 12 '14 at 19:54
see reference.wolfram.com/mathematica/tutorial/… for the options... –  Nasser Jun 12 '14 at 20:01
Thanks, I tried this. There is no more error message but the curve plotted is still weird. It begins oscillating (when e is around -130k) rather than being a smooth curve. Any methods to deal with this? –  riorita Jun 12 '14 at 20:15