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/ > >