Difficulty using NDSolve

Question

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.

Answer 1

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

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