# Stiff system problems with NDSolve

I'm trying to solve the two order differential equations using NDSolve. Here is the code:

ClearAll["Global*"]; Ω1 = 5; Ω2 = -2; m \
= 0.1109; g = 0.0089; timelimit = 10000;
s = NDSolve[{m x''[t] + x'[t] ==
Sin[x[t] - Ω1 t] Exp[-y[t]] +
Sin[2 x[t] - Ω2 t] Exp[-2 y[t]],
m y''[t] + y'[t] ==
Cos[Ω1 t - x[t]] Exp[-y[t]] +
Cos[Ω2 t - 2 x[t]] Exp[-2 y[t]] - g, x[0] == 0,
y[0] == 0, x'[0] == 0, y'[0] == 0}, {x, y}, {t, 0, timelimit}];
ParametricPlot[Evaluate[{x[t], y[t]} /. s], {t, 0, timelimit},
Frame -> True, PlotRange -> All, AspectRatio -> 1/2]


The parameter of m and g is based on physical model. And the source in the right hand of the function is the product of Sin function and exponential function.

If I set y[0] == 0

It will illustrate that NDSolve::ndsz: At t == 3.0092307260023015, step size is effectively zero; singularity or stiff system suspected

If I changed the value of y'[0] to 1,2 or other values, or I add whenevent like this , this problems will be solved.

WhenEvent[y[t] == 0, y'[t] -> - y'[t]]} .

But I don't know the basic mechanism.

LIke in what conditions it will show a stiff system. If I want to get the simulation results in that condition, what's the proper way to solve the problem.

Any suggestions would be appreciated.

• You tagged the question "singularity" but you say "stiff." Which is it, stiffness or a singularity? Related: mathematica.stackexchange.com/questions/39028/… Oct 8, 2021 at 16:08
• Exponential function in the source term is a problem. The system diverges when y(t) becomes negative. Are you sure your equations are physically correct? Oct 8, 2021 at 19:12
• @ Michael E2 I thought singularity may lead to stiff system, so I added it. It may be wrong. Oct 10, 2021 at 16:14
• Yes, true. If I add whenevent and keep y[t]>0, this problem may be sovled Oct 10, 2021 at 16:17

It may help to look at

Plot[{x[t] /. s, y[t] /. s}, {t, 0,  Last[First[x["Domain"] /. First[s]]]}]


It appears that the derivatives are getting very steep. Indeed, the plot

Plot[{x'[t] /. s, y'[t] /. s}, {t, 0,
Last[First[x["Domain"] /. First[s]]]}, PlotLegends -> {"x'", "y'"}]


suggests a singularity at around 3

• Good point. Find singularity through output. What does this {t, 0, Last[First[x["Domain"] /. First[s]]]} mean. Oct 10, 2021 at 16:35
• NDSolve returns InterpolatingFunctions, in this case x and y. They tell you the range over which they are defined. So, for example, x["Domain"]/.First[s] (the first here is because NDSolve returns a list of solutions, we want to look at the first --and I suspect only -- solution in that list) returns the domain {{from_t, to_t}}. Last[First[...]] returns to_t, the maximum t of the domain. Beyond that, the interpolating function relies on extrapolation. Oct 11, 2021 at 18:11
• Oh, I see, thanks @JV3. Recently, I found that sometimes if I increase the total simulation time, it may run into the stiff system. Do you know why the total simulation time will have an effect on the final result? Oct 13, 2021 at 13:23
• Again, I'd recommend looking at the solution to see what's happening at the point where you run into trouble. Oct 14, 2021 at 15:11