[Date Index]
[Thread Index]
[Author Index]
Re: NDSolve for Newtonian Orbits Question
*To*: mathgroup at smc.vnet.net
*Subject*: [mg47304] Re: NDSolve for Newtonian Orbits Question
*From*: "Kevin J. McCann" <kjm at KevinMcCann.com>
*Date*: Mon, 5 Apr 2004 05:22:36 -0400 (EDT)
*References*: <c4jasf$dpk$1@smc.vnet.net>
*Sender*: owner-wri-mathgroup at wolfram.com
David,
You are the one who almost always answers our questions; so, I am pleased to
offer some help. One of the problems with orbit calculations is that pesky
1/r term. NDSolve does not know that r should always be positive. I am not
sure that this is your problem, but I have run into it on several occasions.
As a result, I always do these calculations in Cartesian coordinates. I have
a nb with calculations for Jupiter's Trojan points, if you would like, I can
send it to you.
Kevin
"David Park" <djmp at earthlink.net> wrote in message
news:c4jasf$dpk$1 at smc.vnet.net...
> Dear MathGroup,
>
> I am trying to obtain a numerical solution for Newtonian orbits. (I can
solve the
> du/dphi equation symbolically with DSolve to obtain the equation of an
ellipse, but I also want to know how to do numerical solutions.) I'm having
difficulty in knowing how to use NDSolve. For a start I just want to solve
for the radius r as a function of time.
>
> Here are the equations. We have an effective potential given by
>
> Clear[r]
> Veff[M_, h_][r_] = -M/r + (1/2)*(h^2/r^2)
>
> M is the attracting mass and h is the angular momentum. The following is a
plot of a particular case.
>
> Plot[Veff[1, 1][r], {r, 0.5, 10},
> PlotRange -> All,
> Frame -> True,
> FrameLabel -> {r, Veff},
> Axes -> False,
> ImageSize -> 450];
>
> The differential equation is given by...
>
> deqn[M_, h_, En_] = En == (1/2)*Derivative[1][r][t]^2 + Veff[M, h][r[t]]
> En == h^2/(2*r[t]^2) - M/r[t] + (1/2)*Derivative[1][r][t]^2
>
> where En is the energy. If I pick the energy to be -0.3 and put r[0]
somewhere in the orbit range I should obtain a nice elliptical orbit. But
when I try to solve I run into all kinds of problems.
>
> Clear[r]
> NDSolve[{deqn[1, 1, -0.3], r[0] == 1.5}, r, {t, 0, 10}]
> r[t_] = r[t] /. Part[%, 2]
>
> Plot[r[t], {t, 0, 5.25}];
>
> It appears that the numerical solution gets stuck at the turning points.
How can I obtain an extended solution?
>
> David Park
> djmp at earthlink.net
> http://home.earthlink.net/~djmp/
>
>
Prev by Date:
**Re: Help: Generalized convolution using ListConvolve**
Next by Date:
**Open letter: Who steals memory? :-)**
Previous by thread:
**Re: NDSolve for Newtonian Orbits Question**
Next by thread:
**Re: NDSolve for Newtonian Orbits Question**
| |