How to remove the warnings and errors from this manipulate code?

Question

I'm experiencing several problems with the manipulate code below. This code numerically solves three first-order non-linear differential equations, and then output a 3D plot of the curve, using the initial conditions entered by the user from the three sliders :

Clear["Global`*"]
Theta1 = -20;
Theta2 = 20;
dynamics[om_, or_, ol_] := dynamics[om, or, ol] = NDSolve[{
    oM'[theta] == (oM[theta] + 2oR[theta] - 2oL[theta] - 1)oM[theta],
    oR'[theta] == (oM[theta] + 2oR[theta] - 2oL[theta] - 2)oR[theta],
    oL'[theta] == (oM[theta] + 2oR[theta] - 2oL[theta] + 2)oL[theta],
oM[0] == om, oR[0] == or, oL[0] == ol

}, {oM, oR, oL}, {theta, Theta1, Theta2},
Method -&gt; &quot;StiffnessSwitching&quot;
(* constraints to be imposed : oM, oR, oL &gt; 0 only *)

]
Tmin[om_, or_, ol_] := Theta1 (* to be fixed )
Tmax[om_, or_, ol_] := Theta2 ( to be fixed *)
solution[om_, or_, ol_] := ParametricPlot3D[
    Evaluate[{oM[theta], oR[theta], oL[theta]}/.dynamics[om, or, ol]],
    {theta, Tmin[om, or, ol], Tmax[om, or, ol]}
]
Manipulate[
    Show[
    solution[om, or, ol],
    PlotRange -> {{0, 2}, {0, 2}, {0, 2}},
    SphericalRegion -> True
    ],
    {{om, 0.3, a}, 0, 2, 0.01},
    {{or, 0.0, b}, 0, 2, 0.01},
    {{ol, 0.7, c}, 0, 2, 0.01}
]

At compilation, I'm getting several error messages that I don't know how to solve. The manipulate box should be regular for all values entered from the sliders (from 0 up to 2 or more), especially when any of the variables is set to 0. The oM, oR and oL variables should only be positive, so a constraint should be added to the NDSolve part, and the curve extremities should be properly defined with the Tmin and Tmax definitions above. Currently, it doesn't work well.

How can I fix these problems ?

EDIT : Here are two typical messages that annoys me :

... step size is effectively zero; singularity or stiff system suspected

... lies outside the range of data in the interpolating function. Extrapolation will be used.

You only get errors like InterpolatingFunction::dmval: "Input value {-20.} lies outside the range of data in the interpolating function. Extrapolation will be used. ", right? Or other errors? Always state what the errors are in the OP. — Marius Ladegård Meyer, Dec 27 '16 at 20:50
Also, since this is an initial value problem, why don't you specify the function value at the initial points, e.g. oM[Theta1] == om? — Marius Ladegård Meyer, Dec 27 '16 at 20:51
@MariusLadegårdMeyer, there are several sorts of warning or errors while compilating this code. Yes, I want to prevent any extrapolation. Also, I can't use oM[Theta1] == om since the "initial" conditions are defined specifically at theta = 0, and I need to get the curve on both sides : theta < 0 (backward in time) and theta > 0 (forward in time). — Cham, Dec 27 '16 at 23:56
Related: http://mathematica.stackexchange.com/a/91268/4999 -- Alternatively, stay within the domain. — Michael E2, Dec 28 '16 at 01:29
Please include the message name too (NDSolve::ndsz or whatever they are). It's useful for searching. -- The first one can be taken care of with Quiet[..., NDSolve::ndsz] or whatever, assuming it blows up in finite time, which is not an error. — Michael E2, Dec 28 '16 at 01:31
How can I find the extremities of the curve, without any extrapolation ? How to properly define the end parameters Tmin[om, or, ol] and Tmax[om, or, ol] ? — Cham, Dec 28 '16 at 01:57
See (87525). Also (28337). These, too, are related: (15953), (20774). — Michael E2, Dec 28 '16 at 05:43
@MichaelE2, sorry, these questions/answers aren't helping me. I need to solve this issue : how to properly define the argument (theta min and theta max) of both extremities of the curve, to remove the extrapolation warning ? — Cham, Dec 28 '16 at 16:39
@MichaelE2, well, I don't understand that answer. Please, could you be more specific, and apply your solution as an answer here ? — Cham, Dec 28 '16 at 16:58
I'm sorry I'm leaving on a trip in a few minutes. The key things are "Domain" and "ExtrapolationHandler" in that post. You can also find them in other posts. Maybe someone else will come along in the meantime. — Michael E2, Dec 28 '16 at 17:23
Does Dynamic\Manipulate have anything to do with the question? Can you pose a specific example with specific parameters that demonstrate the error in a more simple way? — george2079, Dec 28 '16 at 18:30
@george2079, already at compilation of the code above, I'm getting several messages, starting with this : Input value [19.9983] lies outside the range of data in the interpolating function. Extrapolation will be used. As soon as I move a bit one of the sliders, I may get more messages, and some like this : At theta == 0.402456, step size is effectively zero; singularity or stiff system suspected. (try these inputs : a = 1.49, b = 1.05 and c = 0.06). — Cham, Dec 28 '16 at 18:44

george2079 · Accepted Answer · 2016-12-28T20:03:46.497

1

a bit of an extended comment, this is a minimal demonstration of the error you see. (You should have posted something like this to begin with) Note dynamic/manipulate has nothing to do with it, nor does the initial condition being in the middle of the domain cause the problem.

r = With[{om = 1, or = 2, ol = .3}, First@NDSolve[{
      oM'[ theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 1) oM[theta],
      oR'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 2) oR[theta],
      oL'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] + 2) oL[theta],
      oM[0] == om, oR[0] == or, oL[0] == ol},
     {oM, oR, oL},
     {theta, 0, 2}]];

NDSolve::ndsz: At theta == 0.35134263731177495`, step size is effectively zero; singularity or stiff system suspected. >>

Plot[{(oR /. r)[x], (oM /. r)[x], (oL /. r)[x]}, {x, 0, 1}, 
 PlotRange -> {0,40}]

the pde is simply blowing up.

increasing ol to .4 yields a stable solution:

Now that we see whats going on, go back to Manpulate, adding Quiet to surpress the error and extract the valid range from each solution:

Manipulate[
 Quiet[
  r = With[{om = 1, or = 2}, First@NDSolve[{
       oM'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 1) oM[
          theta],
       oR'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] - 2) oR[
          theta],
       oL'[theta] == (oM[theta] + 2 oR[theta] - 2 oL[theta] + 2) oL[
          theta],
       oM[0] == om, oR[0] == or, oL[0] == ol},
      {oM, oR, oL},
      {theta, 0, 2}]]];
  range = Flatten[(oM /. r)["Domain"]];
  Plot[{(oR /. r)[x], (oM /. r)[x], (oL /. r)[x]}, 
      Evaluate[Join[{x}, range]], PlotRange -> {0, 40}],
{{ol, .4}, .3, 1}]

Finally, fixing the original should look like this:

solution[om_, or_, ol_] := (
  result = Quiet@First@dynamics[om, or, ol];
  range = (oM /. result)["Domain"][[1]];
  ParametricPlot3D[{oM[theta], oR[theta], oL[theta]} /. result, 
   Evaluate[Join[{theta}, range]]])

edited Dec 28 '16 at 20:03

answered Dec 28 '16 at 19:15

george2079

38,913
1
43
110

I'm not sure it should blow up or not. My 3D graphics is okay, but the numerical resolution should stop if there's really a blow up, and the Tmin, Tmax should take care of it. What do you mean by "increasing ol to .4" ? – Cham Dec 28 '16 at 19:21
see edit, I show how to get the valid range from the solution in case where it blows up. – george2079 Dec 28 '16 at 19:36
With your last code, I'm getting this error message : Join::heads: Heads List and MinMax at positions 1 and 2 are expected to be the same. – Cham Dec 28 '16 at 19:44
1

MinMax is a new function (v10.1), for older versions do: range = (Flatten[(oM /. r)["Grid"]])[[{1, -1}]] .. actually using "Domain" is even better, I updated the code to use that. – george2079 Dec 28 '16 at 19:49
Yep, this last one works very well. Now I'm not sure to understand your code. How do I change my limits Tmin and Tmax in my code above ? – Cham Dec 28 '16 at 19:51
see last edit.. – george2079 Dec 28 '16 at 20:04
It's working great. Sadly, I don't understand that solution ! :-( – Cham Dec 28 '16 at 20:07
While I'm not understanding clearly your last code, I may adopt it anyway since it appears to work flawlesly. I'll wait a bit before marking your answer, since others may have another solution, easier to understand for me. – Cham Dec 28 '16 at 21:12

How to remove the warnings and errors from this manipulate code?

1 Answers1