Legal Notice, Disclaimer, and Book Ordering Information:

(* Legal Notice, Disclaimer, and Book Ordering Information:

Last Modified: FIO January 1, 1995

Mathematica For Physics
Robert L. Zimmerman and Fredrick I. Olness
MathSource Number: 0206-862 

ISBN 0-201-53796-6
For ordering information, call 1-800-822-6339
Addison-Wesley Publishing Company

Communication with the authors:
Fredrick I. Olness:  olness@mail.physics.smu.edu
Robert L. Zimmerman: bob@zim.uoregon.edu

Copyright 1995, Addison-Wesley Publishing Company, Inc. The material in
this file may be distributed freely  so long as the content remains
unchanged, and the copyright and  reference notices are included. The
publisher grants permission for the noncommercial use of these programs
and program segments.  All other uses require the prior written consent
of the publisher.    

All rights are reserved with regard to the  material in the
"Mathematica For Physics" text, and may not be reproduced, stored in a
retrieval system, or transmitted,  in any form or by any means,
electronic, mechanical, photocopying,  recording, or otherwise, without
the prior written permission  of the publisher. 

The programs and applications presented in this book have been included
for  their instructional value. They have been tested with care but are
not  guaranteed for any particular purpose. The publisher does not
offer any  warranties or representations, nor does it accept any
liabilities with  respect to the programs or applications. 

 *)

Chapter 4, Section 2, Problem 4:   
 Numerical Solution for Orbits with Central Forces

<<mainDir:subDir:ch4init.m  (* Mac  Version *)

<<mainDir/subDir/ch4init.m  (* Unix Version *)

Solution

Clear["Global`*"];

eq1= u''[phi]==-u[phi]- m f[1/u[phi]]/(el^2 u[phi]^2);

forceRule= {  V -> (-k/#&)
             ,f[1/u[phi]] -> -V'[r]
             ,r->1/u[phi]
           };

eq2= eq1 //. forceRule

values ={k->1,el->1,m->1 };

initial1={u[0]==1.6,u'[0]==0};

ndsol1= 
NDSolve[{eq2/.values,initial1}//Flatten
          ,u[phi],{phi,0,2 Pi}]//Flatten;

u1[phi_]=u[phi] //.ndsol1;
u1[1.1]

PolarPlot[1/u1[phi] //Evaluate,{phi,0,2 Pi}];

initial2= {u[0]==3,u'[0]==0};

ndsol2= 
NDSolve[{eq2/.values,initial2}//Flatten
          ,u[phi],{phi,0,2 Pi}]//Flatten;

u2[phi_]=u[phi] //.ndsol2;
u2[1.1]

PolarPlot[1/u2[phi] //Evaluate,{phi,0,2 Pi}];

Clear[numOrbit];
numOrbit[V_,r_,r0_:1,dr0_:0,el_:1,m_:1,angle_:(6 Pi)]:=
 Module[{u,phi,uforce,eq1,initial,eqs,ndsol},
   uforce=-D[V,r]/.r->1/u[phi];
   eq1= (u''[phi]== 
            -u[phi]-m uforce/(el^2 u[phi]^2)
          //ExpandAll);
   initial={u[0]==1/r0,u'[0]==-dr0/r0^2};
   eqs = Join[{eq1},initial];
   ndsol= NDSolve[eqs ,u[phi],{phi,0,angle}]//Flatten;  
   PolarPlot[1/u[phi] /.ndsol //Evaluate 
    ,{phi,0,angle}]   
]

numOrbit[-1/r, r, 0.6];

numOrbit[-1/r^1.1,r,0.6];

numOrbit[-1/r^0.9,r,0.6];

numOrbit[-1/r - 0.02/r^3,r,0.6];

numOrbit[-1/r + 0.02/r^3,r,0.6];