0
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Marius Ladegård Meyer Dec 27 '16 at 20:50
  • $\begingroup$ 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? $\endgroup$ – Marius Ladegård Meyer Dec 27 '16 at 20:51
  • $\begingroup$ @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). $\endgroup$ – Cham Dec 27 '16 at 23:56
  • $\begingroup$ Related: mathematica.stackexchange.com/a/91268/4999 -- Alternatively, stay within the domain. $\endgroup$ – Michael E2 Dec 28 '16 at 1:29
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Dec 28 '16 at 1:31
1
$\begingroup$

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

enter image description here

the pde is simply blowing up.

increasing ol to .4 yields a stable solution:

enter image description here

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]]])
$\endgroup$
  • $\begingroup$ 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" ? $\endgroup$ – Cham Dec 28 '16 at 19:21
  • $\begingroup$ see edit, I show how to get the valid range from the solution in case where it blows up. $\endgroup$ – george2079 Dec 28 '16 at 19:36
  • $\begingroup$ 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. $\endgroup$ – Cham Dec 28 '16 at 19:44
  • 1
    $\begingroup$ 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. $\endgroup$ – george2079 Dec 28 '16 at 19:49
  • $\begingroup$ 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 ? $\endgroup$ – Cham Dec 28 '16 at 19:51

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.