Stop Integration in NDSolve

The following code executes properly for some choices of initial points and parameter c, but fails for other choices. Could the failure be due to computation outside of the defined t, x, and y ranges? If so, I don’t know how to stop the computation with a WhenEvent option.

Manipulate[
tmin = -20; tmax = 20;
xmin = -3 Pi; xmax = 3 Pi;
ymin = -4.5; ymax = 4.5;

z[t_, c_, x0_, y0_] := {x[t], y[t]} /.
First@NDSolve[{x'[t] == y[t], y'[t] == -Sin[x[t]] – c*y[t],
x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, tmin, tmax}];

ClickPane[
Show[ ParametricPlot[g, {t, tmin, tmax},
PlotRange -> {{-3.2 Pi, 3.2 Pi}, {-4.3, 4.3}}, Frame -> True,
Axes -> None, ImageSize -> 700, PlotStyle -> Black],
vf, Epilog -> {PointSize[Large], Point[sp]}],

(AppendTo[g, z[t, c, #[[1]], #[[2]]]]; {t0, x0} = #;
AppendTo[sp, #]) &],

{{c, 0.5, “c”}, 0, Sqrt[2], 0.01},

Row[{Button[“Delete all integral curves”, {g = {{}, Axes -> None};
sp = {}},ImageSize -> {150, 20}], Spacer[20],
Button[“Reset c”, {c = 0.5}, ImageSize -> {50, 20}]}],

SaveDefinitions -> True,

Initialization :> {(g = {{}, Axes -> None};
sp = {}; {x0, t0} = {-1, -1}),

vf := VectorPlot[{y, -Sin[x] – c*y}, {x, xmin, xmax}, {y, ymin,
ymax}, VectorPoints -> {22, 15}, PerformanceGoal -> “Quality”,
VectorScale -> {0.03, 1.3, None}, PerformanceGoal -> “Speed”]
}
]

With c=0.5 click on a point in the window to get a solution through that point. After about 6 clicks, the solution is not plotted. Two error messages are posted:
1. NDSolve::dsvar: -19.9992 cannot be used as a variable.
2. ReplaceAll::reps: (yadayadayada….. )

=================

  

 

Simply move the delayed definition of z outside of Manipulate, and everything will work just fine. Also, you will notice that this will stop the continuous updating that was plaguing your Manipulate before. It will only update when a mouse click is received.
– MarcoB
Jun 15 at 21:39

=================

2 Answers
2

=================

As I mentioned in comments, your main problem was the placement of the definition of z; that definition does not need to be re-evaluated every time there is a modification within the body of Manipulate since it’s already delayed. The correct placement for such a definition would be in the Initialization code for the Manipulate. Similarly, the definitions of your constants only need to be evaluated once, so I moved them to the initialization portion as well.

Here is also a somewhat cleaned up version of your code; your original code seemed to contain a lot of leftover assignments and conflicting options that did not have any effect, perhaps left over from another application:

Manipulate[
ClickPane[
Show[
ParametricPlot[
g, {t, tmin, tmax}, PlotRange -> {{-3.2 Pi, 3.2 Pi}, {-4.3, 4.3}},
Frame -> True, Axes -> None, ImageSize -> 700, PlotStyle -> Black
],
vf,
Epilog -> {PointSize[Large], Point[sp]}
],
(AppendTo[g, z[t, c, Sequence @@ #]]; AppendTo[sp, #]) &
],

(* Manipulate variable *)
{{c, 0.5}, 0, Sqrt[2], 0.01},

(* reset buttons *)
Row[{
Button[“Delete all integral curves”, g = {}; sp = {}, ImageSize -> {150, 20}],
Spacer[20],
Button[“Reset c”, c = 0.5, ImageSize -> {50, 20}]
}],

(* One-time initialization code *)
Initialization :> (
g = {};
sp = {};
vf := VectorPlot[
{y, -Sin[x] – c*y}, {x, xmin, xmax}, {y, ymin, ymax},
VectorPoints -> {22, 15}, PerformanceGoal -> “Quality”,
VectorScale -> {0.03, 1.3, None}
];
z[t_, c_, x0_, y0_] := {x[t], y[t]} /. First@
NDSolve[
{x'[t] == y[t], y'[t] == -Sin[x[t]] – c*y[t], x[0] == x0, y[0] == y0},
{x[t], y[t]}, {t, tmin, tmax}
];
tmin = -20; tmax = 20;
xmin = -3 Pi; xmax = 3 Pi;
ymin = -4.5; ymax = 4.5;
)
]

With these changes, everything works well:

  

 

@ MarcoB Thanks for the cleaned up code; it has fixed the problem that I thought was generated by the variables out of range. But the other problem persists. Set c =0.8 or even larger and click away from the x-axis. If you click very close to the x-axis, an integral curve will be plotted. If you click too far away, only the initial point will be plotted. This is strange. Can it be an issue with NDSolve?
– Stephen
Jun 15 at 22:15

  

 

I didn’t see your latest code before I sent off my comment a few moments ago. BTW, R U sure that that your posted code is correct? There may be a typo as I cant get it to work
– Stephen
Jun 15 at 22:21

  

 

@Stephen Thank you for pointing that out; I was missing a semicolon before the definition of z in the initialization, due to careless copy/paste. It’s fixed now. As for the high c issue, have you tried to calculate a solution outside of manipulate with that value of c, to see if a solution is returned by NDSolve in those conditions?
– MarcoB
Jun 15 at 22:44

  

 

MarcoB: The non-plotting issue has to do with the t-interval. If I shorten it sufficiently, the plot will show. Since this ODE models a nonlinear damped system (e.g., a pendulum), when c (the damping coefficient gets large), the system not only goes to zero very rapidly as t->+infinity, but it also goes becomes unbounded rapidly as t->-infinity. I can finesse the issue by requiring tmin=0 or by implementing WhenEvent[Abs[y[t]] > 4,”StopIntegration”]
– Stephen
Jun 16 at 3:12

  

 

I changed the time interval to be -5<=t<=100. That takes care of all issues. When t becomes too negative,y is out of range of the extrapolating function – Stephen Jun 16 at 3:30 I have revised the code to use StreamPlot instead of VectorPlot and NDSolve. The latter command created issues when the domain for the solution extended to negative times. I don't know why, but NDSolve produced integral curves that did not follow the direction field produced by VectorPLot as provided by MarcoB's nice reformulation of my original code. I will continue to work on this issue. In the meantime, StreamPlot provides an easy fix. Manipulate[ Show[StreamPlot[{v , -Sin[u] - c v}, {u, -3 Pi, 3 Pi}, {v, -4.5, 4.5}, StreamScale -> Automatic,
AspectRatio -> 0.5,
ImageSize -> 700,
Frame -> True,
StreamPoints -> Fine,
Epilog -> {{Red, PointSize[0.01],
Point[{{0, 0}, {-Pi, 0}, {Pi, 0}, {-2 Pi, 0}, {2 Pi, 0}}]}}]
],
{{c, 0.5, “c”}, 0, 3, 0.01},
Button[Style[“Reset”, FontFamily -> “Helv”, 14], {c = 0.5},
ImageSize -> {50, 25}]
]

  

 

Thanks for the follow up! (+1)
– MarcoB
Jun 17 at 5:23