Strange plot results from NDSolve with ExplicitEuler

Question

The following code attempts to plot Euler approximations to an autonomous system of ODEs of the form $x'=y$, $y'=-x$. Click at a point $(x_0,y_0)$ in the $xy$-window to generate a plot based on the ExactEuler method using NDSolve on the time interval [0,$t_f$] with step size $h$ and initial values $x(0)=x_0, y(0)=y_0$

Panel@DynamicModule[{g = {}, p = {}, sol, x0, y0, tf, h},
  sol[{x0_, y0_}, h_] := {x[t], y[t]} /. 
    First@NDSolve[{x'[t] == y[t], y'[t] == -x[t], x[0] == x0, 
       y[0] == y0}, {x, y}, {t, -10, 10}, 
       Method -> {"FixedStep", Method -> "ExplicitEuler"}, 
       StartingStepSize -> h, MaxStepFraction -> 1];
  Column[{
    Row[{Control[{{tf, Pi, Style["tf ="]}, ImageSize -> 40}], Spacer[30], 
    Control[{{h, 1.0, Style["h ="]}, ImageSize -> 40}], Spacer[220], 
    Button[Style["Delete all solutions"], g = {}; p = {}, 
    ImageSize -> {Automatic, 25}]}],
ClickPane[
    Dynamic@Show[
    ParametricPlot[g, {t, -tf, tf}, PlotRange -> {{-2, 2}, {-2, 2}},
    PlotPoints -> Ceiling[tf/h], MaxRecursion -> 0, 
    ImageSize -> 400, Axes -> None, Frame -> True, 
    PlotLabel -> Dynamic[Style[MousePosition["Graphics"]]]],
    Graphics[{PointSize[Large], Point[p]}]],
   (AppendTo[g, sol[#, h]]; AppendTo[p, #]) &]}]]

The strange plots are illustrated below. Start with the upper left plot; the solution through the point $(0,0.5)$ with $h=1.0$ on the time interval $[0,2\pi]$ looks reasonable. For the upper right plot set $h=0.5$ and click on the point $(0,0.5)$ to obtain a second solution through $(0,0.5)$. Note how the first solution (with $h=1.0)$ has shifted. The lower right plot shows a third solution through $(0,0.5)$ with $h=0.25$. Note how the first two solutions morphed into new curves. The situation becomes even weirder in the lower right plot when the solution through $(0,0.5)$ with $h=0.1$ as added. Note how the polygonal form of the $h=1.0)$ solution has morphed into a smooth curve.

It gets weirder still when the plotting interval is expanded to go backward in time (by modifying the time interval for ParametricPlot to be $[-t_f,t_f]$. The next image shows plots of "solutions" on $[-\pi,\pi]$ with $h=1.0$ produced by random clicks in the $xy$-plane. Notice how each solution is disjoint from it's initial point.

I suspect that the strange behavior may be result from my improper use of the PlotPoints and MaxRecursion options for ParametricPlot. Although I can code the Euler method directly from its simple formula, I need to stick with the NDSolve method and the same AppendTo construction you see in my code.

I don't know how to resolve these issues.

Answer 1

To plot just the steps, extract the "ValuesOnGrid" from the interpolating functions and plot them with ListLinePlot:

Panel@DynamicModule[{g = {}, p = {}, sol, x0, y0, tf, h}, 
  sol[{x0_, y0_}, h_] := {x[t], y[t]} /. 
    First@NDSolve[{x'[t] == y[t], y'[t] == -x[t], x[0] == x0, 
       y[0] == y0}, {x, y}, {t, -10, 10},
      Method -> {"FixedStep", Method -> "ExplicitEuler"}, 
      StartingStepSize -> h, MaxStepFraction -> 1];
  Column[{
    Row[{
      Control[{{tf, Pi, Style["tf ="]}, ImageSize -> 40}],
      Spacer[30],
      Control[{{h, 1.0, Style["h ="]}, ImageSize -> 40}],
      Spacer[220],
      Button[
       Style["Delete all solutions"],
       g = {}; p = {},
       ImageSize -> {Automatic, 25}]}],
    ClickPane[
     Dynamic@Show[
       ListLinePlot[                                 (* replaces ParametricPlot *)
        Transpose[# /. t -> "ValuesOnGrid"] & /@ g,  (* constructs the points *)
        ImageSize -> 400, Axes -> None, Frame -> True, 
        AspectRatio -> Automatic, 
        PlotLabel -> Dynamic[Style[MousePosition["Graphics"]]]],
       Graphics[{PointSize[Large], Point[p]}],
       PlotRange -> {{-2, 2}, {-2, 2}}               (* move PlotRange here *)
       ],
     (AppendTo[g, sol[#, h]]; AppendTo[p, #]) &]}]]

As J.M. observed, the interpolating functions automatically use cubic Hermite interpolation. The adjustment of the number of points plotted with PlotPoints -> Ceiling[tf/h] means the number of points plotted increases as h decreases, if you follow the protocol in the OP. As the number of plot points becomes much greater than the number of steps, the cubic nature of the interpolation becomes more apparent.