0
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ You tagged the question "singularity" but you say "stiff." Which is it, stiffness or a singularity? Related: mathematica.stackexchange.com/questions/39028/… $\endgroup$
    – Michael E2
    Commented Oct 8, 2021 at 16:08
  • 2
    $\begingroup$ 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? $\endgroup$
    – yarchik
    Commented Oct 8, 2021 at 19:12
  • $\begingroup$ @ Michael E2 I thought singularity may lead to stiff system, so I added it. It may be wrong. $\endgroup$
    – Yue Yu
    Commented Oct 10, 2021 at 16:14
  • $\begingroup$ Yes, true. If I add whenevent and keep y[t]>0, this problem may be sovled $\endgroup$
    – Yue Yu
    Commented Oct 10, 2021 at 16:17

1 Answer 1

Reset to default
2
$\begingroup$

It may help to look at

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

x[t] in blue, y[y] in orange

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'"}]

x'[t] and y'[t]

suggests a singularity at around 3

Improve this answer
$\endgroup$
4
  • $\begingroup$ Good point. Find singularity through output. What does this {t, 0, Last[First[x["Domain"] /. First[s]]]} mean. $\endgroup$
    – Yue Yu
    Commented Oct 10, 2021 at 16:35
  • $\begingroup$ 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. $\endgroup$
    – user46831
    Commented Oct 11, 2021 at 18:11
  • $\begingroup$ 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? $\endgroup$
    – Yue Yu
    Commented Oct 13, 2021 at 13:23
  • $\begingroup$ Again, I'd recommend looking at the solution to see what's happening at the point where you run into trouble. $\endgroup$
    – user46831
    Commented Oct 14, 2021 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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