# 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 == om, oR == or, oL == 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. 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? 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. 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. 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 == om, oR == or, oL == 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 == om, oR == or, oL == 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"][];
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. 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. 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