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

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 -> "StiffnessSwitching"
(* constraints to be imposed : oM, oR, oL > 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: mathematica.stackexchange.com/a/91268/4999 -- Alternatively, stay within the domain. – Michael E2 Dec 28 '16 at 1: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 1:31

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]]])

• 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
• 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