0
$\begingroup$

I am tying to solve two differential equations with two delays numerically. The result didn't come out and it was written that “NDSolve::ndsz:...step size is effectively zero; singularity or stiff system suspected”

Here's the code.

Manipulate[
Module[{plt1, plt2, sol, x0 = xx0, y0 = yy0}, 
sol = NDSolve[{x'[t] == d*y[t] + k*y[t - τ2], 
 y'[t] == μ*y[t - τ2]*(b*(x[t] + β)^2 + c) - μ*
    h*x[t] - μ*a*x[t - τ1] + μ*(f - g), 
 x[t /; t <= 0] == x0, y[t /; t <= 0] == y0}, {x[t], y[t]}, {t, 0,
  50}]; plt1 = ParametricPlot3D[{x[t], y[t]} /. sol, {t, 0, 50}, 
PlotRange -> All, AspectRatio -> 0.6, 
PlotStyle -> {RGBColor[1, 0, 1], Thick}, AxesLabel -> {x, y}];
plt2 = StreamPlot[{d*y + 
  k*y, μ*y*(b*(x + β)^2 + c) - μ*h*x - μ*a*
   x + μ*(f - g)}, {x, 0, 50}, {y, 0, 50}, FrameLabel -> {
  {Style[Row[{Style["x", Italic], "(", Style["t", Italic], ")"}], 
    14], None},
  {Style[Row[{Style["y", Italic], "(", Style["t", Italic], ")"}], 
    14], Style[Row[{Style["The phase portrait of system", Bold]}],
     14]}}, StreamPoints -> 50];
start = Graphics3D[Point3DBox[{xx0, yy0}]];
plt3 = {plt1, start};
Show[plt3, ImageSize -> {450, 400}]],
Style["Persamaan differensial :", Bold], Style["\!\(\*OverscriptBox[\(\(x\)\(\\\ \)\), \(.\)]\)= \d*y[t]+k*y[t-τ2] ", Bold],
Style["\!\(\*OverscriptBox[\(y\), \(.\)]\) = \μ*y[t-τ2]*(b*(x[t]+β\!\(\*SuperscriptBox[\()\), \(2\)]\\)+c)-μ*h*x[t]-μ*a*x[t-τ1]+μ*(f-g) ", Bold], 
Delimiter,
Style["parameters", Bold, 10],{{d, 1, "d"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{β, -2, "β"}, -100, 0, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{c, -2, "c"}, -100, 0, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{μ, 0.13, "μ"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{a, 12, "a"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{g, 20, "g"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{h, 10, "h"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{f, 60, "f"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{b, 1, "b"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{k, 1, "k"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{τ1, 0, "τ1"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"},
{{τ2, 1, "τ2"}, 0, 100, .01, ImageSize -> Small, Appearance -> "Labeled"}, Delimiter,
Style["initial conditions", Bold,10], {{xx0, 0.4, "\!\(\*SubscriptBox[\(x\), \(0\)]\)"}, 0, 50, .01, ImageSize -> Small,Appearance -> "Labeled"}, 
{{yy0, 0.4,"\!\(\*SubscriptBox[\(y\), \(0\)]\)"}, 0, 50, .01,
ImageSize -> Small, Appearance -> "Labeled"}, ControlPlacement -> Left, SynchronousUpdating -> False]`

I have tried the methods that have been posted in the web but it didn't work. Please give me answer how to deal with it. Thanks in advance for your help.

$\endgroup$
  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful $\endgroup$ – Michael E2 Mar 13 '16 at 12:58
  • 1
    $\begingroup$ Does your DE have a singularity? Method -> "StiffnessSwitching" can handle some forms of stiffness. Search for StiffnessSwitching in the docs. $\endgroup$ – Michael E2 Mar 13 '16 at 13:00
1
$\begingroup$

As noted in the question, the fundamental problem here is with NDSolve. However, the numerous syntax errors in the code make it difficult to analyze this problem as presented. So, extract NDSolve from the code and run it with the default parameters from Manipulate.

sol = NDSolve[{x'[t] == d*y[t] + k*y[t - τ2], 
    y'[t] == μ*y[t - τ2]*(b*(x[t] + β)^2 + c) - μ*h*x[t] - μ*a*x[t - τ1] + μ*(f - g), 
    x[t /; t <= 0] == x0, y[t /; t <= 0] == y0} /. {μ -> 0.13, d -> 1, k -> 1, τ1 -> 0, 
    τ2 -> 1, b -> 1, c -> -2, β -> -2, a -> 12, f -> 60, g -> 20, h -> 10, 
    x0 -> 0.4, y0 -> 0.4}, {x[t], y[t]}, {t, 0, 50}]

It fails at t = 5.07 with the error message given in the question. Using Method -> "StiffnessSwitching", as suggested by Michael E2, often works in such cases but not here. These nonlinear, delay differential equations may be intrinsically singular at finite t for the parameters chosen. However, a slight change to the parameters, namely τ2 -> .067, does produce a stable result out to t = 50.

tm = ((x[t] /. sol) /. t -> "Domain")[[1, 1, 2]]
(* 50. *)
Plot[{x[t] /. sol, y[t] /. sol}, {t, 0, tm}, AxesLabel -> {t, "x,y"}]

enter image description here

Note that the upper bound of the plot, here 50 is extracted from the NDSolve solution, so that Plot can provide a reasonable result, even if NDSolve terminates sooner than t = 50. For instance, the identical code with the original τ2 -> 1 yields

enter image description here

$\endgroup$
  • $\begingroup$ @Aulia, If you plot the default solution (ListLinePlot[{x, y} /. First@sol] is easy, if the OP changes the call to {x, y} instead of {x[t], y[t]}), the plot shows a probable singularity. You can experiment whether this is a numerics problem or a DE problem by setting WorkingPrecision ->16 and higher. If the failure is at 5.07, then it's probably a real singularity. You can write code to deal with the singularity (look up Check and Quiet for instance); a singularity is not necessarily a bad thing programatically. $\endgroup$ – Michael E2 Mar 13 '16 at 15:54
  • 2
    $\begingroup$ I got an error I have never seen: If I set τ2 -> .067 and also change the interval of integration to {t, 0, Infinity}, I get NDSolve::ndssc: Step size changed sign at.... (Nothing to do with your answer, but a new error seemed fun to me. Fixed with MaxStepSize. I was doing to it to find where the solution blows up. Just having fun.) (+1) $\endgroup$ – Michael E2 Mar 13 '16 at 16:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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