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.
NDSolve[]
is performing cubic Hermite interpolation by default. Unfortunately, even when I tried explicitly settingInterpolationOrder -> 1
inNDSolve[]
, I still got a cubic interpolant. Hmm... $\endgroup$NDSolve
. $\endgroup$