Re: submission of a new problem; NDSolve artifacts; version 6.0.1
- To: mathgroup at smc.vnet.net
- Subject: [mg87021] Re: submission of a new problem; NDSolve artifacts; version 6.0.1
- From: "Szabolcs HorvÃt" <szhorvat at gmail.com>
- Date: Sat, 29 Mar 2008 04:23:28 -0500 (EST)
- References: <fsib5n$5ug$1@smc.vnet.net> <47ECFF4D.1000301@gmail.com>
No, you misunderstood me. Please do not send the notebook only to me.
Instead, strip out *everything* from the commands that is not
essential for reproducing the problem. Then send the set of commands
to MathGroup too, as plain text. (I am CCing this message to
MathGroup.) Do make sure that the reduced set of commands can be
evaluated on its own, and does not depend on any variables/function
that you did not provide.
You have to show that you have tried to investigate the problem, and
at least removed the non-essential details. Please do this next time.
I was expecting something like this:
tf = 45.3
tstart = 0;
\[Mu] = 1/2;
vxinitial = -.029656755023561942249954448413973295;
solution[orbtime_, Vx0_, Vy0_, \[Mu]_] :=
NDSolve[
{x''[t] ==
2 y'[t] +
x[t] - ((1 - \[Mu]) (x[t] + \[Mu]))/((x[t] + \[Mu])^2 + y[t]^2)^(
3/2) - (\[Mu] (x[t] - 1 + \[Mu]))/((x[t] - 1 + \[Mu])^2 +
y[t]^2)^(3/2),
y''[t] == -2 x'[t] +
y[t] - ((1 - \[Mu]) y[t])/((x[t] + \[Mu])^2 + y[t]^2)^(
3/2) - (\[Mu] y[t])/((x[t] - 1 + \[Mu])^2 + y[t]^2)^(3/2),
x[0] == 0, y[0] == 858/1000, x'[0] == Vx0, y'[0] == Vy0},
{x[t], y[t]},
{t, 0, orbtime},
Method -> "StiffnessSwitching"]
ParametricPlot3D[
Evaluate[{x[t], y[t], 0} /. solution[tf, vxinitial, 0, 1/2]], {t,
tstart, tf}]
Now the solution:
If I understand correctly, the problem is the strange line that
appears on the plot. This artefact is not from NDSolve, but from
ParametricPlot3D, whose adaptive sampling method interacts badly with
the interpolating function returned from NDSolve. It skips quite a
large range in t. To see this, get the t values at which
ParametricPlot3D samples the function and create a plot manually:
sol = First@solution[tf, vxinitial, 0, 1/2];
pts = Reap[
ParametricPlot3D[
Evaluate[{x[t], y[t], 0} /. sol], {t, tstart, tf},
EvaluationMonitor :> Sow[t]]][[2, 1]];
Manipulate[
ListPlot[Take[Transpose[{x[t], y[t]} /. sol /. t -> Sort[pts]], n],
Joined -> True],
{n, 1, Length[pts]}
]
As the slider is moved, the line "jumps" at one point (more precisely
between the 444th and 445th points):
With[{diff = Differences[Sort[pts]]},
Position[diff, Max[diff]]]
This has been a common problem even with Plot in earlier versions of
Mathematica, but I rarely see it in version 6. It can be cured by
playing with some of the options of ParametricPlot3D, such as
PlotPoints and MaxRecursion.
For example, the following works:
ParametricPlot3D[Evaluate[{x[t], y[t], 0} /. sol], {t, tstart, tf},
PlotPoints -> 30]
On Fri, Mar 28, 2008 at 3:52 PM, Robert M. Lurie <RMLURIE at alum.mit.edu> wrote:
> I have attached the actual notebook that I copied in the previous email.
> There is not a single Mathematica command that causes the problem, but
> the actual conditions in the NDSolve operations.
>
> The final time tf will show the artifact for certain values ( such as
> 45.3) but not for other values (such as 46.0).
> I have added to this notebook the plots of x and y which indicate some
> sharp reversals of direction -- is this a potential answer?
> I welcome any suggestions.
> Many thanks,
> Robert M Lurie
>
>
> Szabolcs Horvát wrote:
> > Robert M. Lurie wrote:
> >> I have been using NDSolve as used in
> >> "Restricted Three-Body Problem in 3D" from The Wolfram
> >> Demonstrations Project
> >> http://demonstrations.wolfram.com/RestrictedThreeBodyProblemIn3D/
> >> Using Version 6.0.1, I have modified the demonstration to show the
> >> classical planar, equal large masses (mu=1/2) cases as in Gass
> >> "Mathematica for Scientists and Engineers: Using Mathematica to do
> >> Science"
> >> Also see
> >> http://scienceworld.wolfram.com/physics/RestrictedThree-BodyProblem.html
> >> Depending on the values, there will be lines emanating from points
> >> where orbits cross. For example at tf (final time) of 45.3, the
> >> figure below is shown. However at tf = 46 no such artifact is seen.
> >> How can NDSolve be modified to eliminate these artifacts?
> >>
> >
> >
> > It is not possible to give an answer unless you tell us *exactly* how
> > to reproduce the problem. Just post a single Mathematica command that
> > reproduces the problem.
> >
> >
>
>