Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: NDSolve for Newtonian Orbits Question

  • To: mathgroup at smc.vnet.net
  • Subject: [mg47308] Re: NDSolve for Newtonian Orbits Question
  • From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
  • Date: Mon, 5 Apr 2004 05:22:46 -0400 (EDT)
  • Organization: Universitaet Leipzig
  • References: <c4jasf$dpk$1@smc.vnet.net>
  • Reply-to: kuska at informatik.uni-leipzig.de
  • Sender: owner-wri-mathgroup at wolfram.com

Hi,

use the Newtonian equations and not the energy conservation.

Regards
  Jens

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/



  • Prev by Date: combining two contour plots
  • Next by Date: Re: Help: Generalized convolution using ListConvolve
  • Previous by thread: Re: NDSolve for Newtonian Orbits Question
  • Next by thread: Re: NDSolve for Newtonian Orbits Question