4
$\begingroup$

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.

enter image description here

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.

enter image description here

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ The problem seems to be due to the fact that NDSolve[] is performing cubic Hermite interpolation by default. Unfortunately, even when I tried explicitly setting InterpolationOrder -> 1 in NDSolve[], I still got a cubic interpolant. Hmm... $\endgroup$ Commented Aug 22, 2017 at 11:53
  • $\begingroup$ Your second problem (the lines not going through the initial point) probably has to do with the fact that the you have evaluated each curve at an even number of points (4), but the original grid from your Euler integration used an odd number of points. This may be one of those rare times when it is useful to reach under the hood and grab the data underlying the interpolating function that NDSolve returns. $\endgroup$ Commented Aug 22, 2017 at 11:54
  • 1
    $\begingroup$ @J.M. "For an ODE of order $n$,...$2n+1$seems like a lower bound on the interpolation order" of solution obtained from NDSolve. $\endgroup$ Commented Aug 22, 2017 at 12:20

1 Answer 1

Reset to default
1
$\begingroup$

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, #]) &]}]]

Mathematica graphics

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Your answer does not allow me to set the number of steps to use; namely Ceiling[tf/h]. The plot through any point on the time interval [0,2Pi] should approximate a circle. Perhaps I don't understand the role of ValuesOnGrid. $\endgroup$ Commented Aug 22, 2017 at 13:12
  • $\begingroup$ Michael 2E: The four plots were based on the interval {t,0,tf} and not {t,-tf,tf} as posted in my code. Sorry for any confusion this might have caused. Modifying your code to replace {t,-10,10} in NDSolve with {t,0,tf} appears make ListLinePlot generate twice as many steps. $\endgroup$ Commented Aug 22, 2017 at 14:49
  • 1
    $\begingroup$ @Stephen PlotPoints -> Ceiling[tf/h] does not set the number of steps in the NDSolve solution. It sets the initial number of points used by ParametricPlot. (I don't think the points even correspond to the steps in the solutions.) Is that what you want? -- Oops, I overlooked tf (somehow I had thought it was the same as 10). ListLinePlot just plots data, and the data used in my answer is the whole solution, not just the solution up to tf. Sorry about that. I might just delete this answer. I'm not sure what you want to do or why, but I can see this does not do it. $\endgroup$ Commented Aug 22, 2017 at 20:27
  • $\begingroup$ What I want is just to plot an Euler solution from t=0 to t=tf. More generally, how do I extract the appropriate InterpolatingPoints for ListLinePlot? $\endgroup$ Commented Aug 22, 2017 at 20:34
  • 1
    $\begingroup$ @Stephen, in that case, to get a picture of the piecewise linear solution, just use the "ValuesOnGrid" property as Michael does and feed that to ListLinePlot[]; no need for ParametricPlot[] anymore. If you want to be able to interpolate linearly as well, then you will have to use Interpolation[] with InterpolationOrder -> 1 on the values obtained from "ValuesOnGrid". $\endgroup$ Commented Aug 22, 2017 at 20:42

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.