Re: NDSolve for Newtonian Orbits Question
- To: mathgroup at smc.vnet.net
- Subject: [mg47317] Re: [mg47297] NDSolve for Newtonian Orbits Question
- From: Selwyn Hollis <sh2.7183 at misspelled.erthlink.net>
- Date: Mon, 5 Apr 2004 05:22:55 -0400 (EDT)
- References: <200404020831.DAA13114@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi David,
You're running into problems because your equation has (r'[t])^2 taking
on negative values. You can get around the problem by differentiating
each side to obtain a second-order equation. Your value of En would
then enter in via an initial condition. Try this:
Veff[M_, h_][r_] = -M/r + (1/2)*(h^2/r^2)
deqn[M_, h_] = 0 == Simplify[(1/r'[t])*D[(1/2)*r'[t]^2 + Veff[M,
h][r[t]], t]]
NDSolve[{deqn[1,1], r[0]==1.5, r'[0]==Sqrt[2*(-0.3-Veff[1, 1][1.5])]},
r, {t,0,10}]
rr[t_] = r[t] /. First[%]
Plot[rr[t], {t, 0, 10}];
-----
Selwyn Hollis
http://www.math.armstrong.edu/faculty/hollis
(edit reply-to to reply)
On Apr 2, 2004, at 3:31 AM, David Park wrote:
> 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/
>
>
>
- References:
- NDSolve for Newtonian Orbits Question
- From: "David Park" <djmp@earthlink.net>
- NDSolve for Newtonian Orbits Question