{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 1 16 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 16 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 257 1 {CSTYLE "" -1 -1 "Times" 1 16 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT -1 14 "Kepler Problem" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 107 "We demonstrate the numerical solution of Newton's equation for the Kepler problem in Cartesian coordinates." }}{PARA 0 "" 0 "" {TEXT -1 192 "A unit system is chosen to make the equations dimensionless, and we assume that the Sun is infinitely heavy, i.e., we just consider the motion of a singl e body around the center of attraction." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 41 "This worksheet contains several sect ions:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 147 "1) In the first section we introduce a system of coordinates that ren der the solar system dimensionless. This is achieved by choosing three scales:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 58 "the AU (average Sun-Earth distance) as a distance measure;" }} {PARA 0 "" 0 "" {TEXT -1 37 "the solar mass M_s as a mass measure;" }} {PARA 0 "" 0 "" {TEXT -1 74 "the year (time it takes for Earth to go a round the Sun) as a time measure." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 167 "We solve the Kepler problem numerically \+ for this problem, and demonstrate a few features, such as the virial t heorem (relation between kinetic and potential energies)." }}{PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 132 "2) We then explor e the problem of radial motion, including a discussion of stability of circular orbits as a function the force law;" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 70 "3) We then conclude with \+ a brief discussion of the scattering problem." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "We begin with sec tion 1:" }}{PARA 0 "" 0 "" {TEXT -1 54 "For the solar system we have t he following parameters:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 24 "Solar mass: 2*10^\{30\} kg" }}{PARA 0 "" 0 "" {TEXT -1 53 "Gravitational constant: 6.67*10^\{-11\} SI (N*m^2/kg^2)" }}{PARA 0 "" 0 "" {TEXT -1 24 "Earth mass: 6*10^\{24\} kg" }}{PARA 0 " " 0 "" {TEXT -1 34 "Sun-Earth distance 1.5 * 10^\{11\} m" }}{PARA 0 " " 0 "" {TEXT -1 32 "Earth's orbit time: 3.15*10^\{7\}s" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 17 "G:=6.67*10^(-11);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "M_s:=1.99*10^(30);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "M_e:=5.97*10^(24);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "G*M_s*M_e;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "AU:=1.5*10^(11);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "G*M _s*M_e/AU^2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 126 "Earth's attracti on towards the Sun is rather large when expressed in Newtons, but the \+ acceleration in m/s^2 is not very large:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "G*M_s/AU^2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 234 " Note that this gravitational acceleration is balanced at the center of the Earth by the centrifugal force due to the rotation about the Sun. The only net effect on an object at the surface of the Earth is in th e form of a tidal force." }}{PARA 0 "" 0 "" {TEXT -1 23 "The period is given as:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "T_e:=evalf(36 5.26*24*60*60);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "v_e:=eva lf(2*Pi*AU)/T_e;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 46 "Earth moves w ith about 30 km/s around the Sun." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 278 "Now we introduce a unit system in which \+ the Earth's orbit radius is of order unity, the masses are measured in multiples of Sun's mass, and times are measured in terms of orbit tim es. We can pick three dimensionful quantities and set the scales as fo llows [remember: SI=MKS(A)]:" }}{PARA 0 "" 0 "" {TEXT -1 12 "M = M' * \+ M_s" }}{PARA 0 "" 0 "" {TEXT -1 12 "T = T' * T_e" }}{PARA 0 "" 0 "" {TEXT -1 11 "X = X' * AU" }}{PARA 0 "" 0 "" {TEXT -1 78 "and with the \+ substitution we express equations in terms of the new variables. " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 419 "A simila r procedure is used in modern physics to scale the atomic orbital prob lem; one notable exception is that one uses the electron mass as a sca le. In the solar Kepler problem this is not done, because the planet's mass drops out of the equation of motion. In the quantum problem it d oes not, because the force of attraction is electric, i.e., not propor tional to the lighter particle's (i.e., the electron's) mass." }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 96 "What is o ur unit of velocity in these natural units for the Earth-around-the-Su n motion problem?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "unitLe ngth:=AU; # introduced unit #1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "unitTime:=T_e; # introduced unit #2" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "unitVel:=unitLength/unitTime; # a derived unit!" } }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 134 "The velocities V' will come ou t as multiples of about 4.7 km/s, which means that Earth's cruising sp eed [about the Sun] becomes 2 Pi!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "v_e/unitVel;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 26 "What is the unit of force?" }}{PARA 0 "" 0 "" {TEXT -1 73 "The natural unit should be the force experienced by \+ Earth due to gravity." }}{PARA 0 "" 0 "" {TEXT -1 78 "The unit for acc eleration is clearly the velocity unit divided by a time unit:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "unitAcc:=unitVel/unitTime; # a derived unit!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "unitMas s:=M_s; # introduced unit #3 !" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "unitForce:=unitMass*unitAcc; # a derived unit!" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 170 "Given this unit for the acceleration, it is of interest to express the acceleration that Earth experiences fro m the gravitational attraction towards the Sun in this unit:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "(G*M_s/AU^2)/unitAcc;" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "The force experienced by Earth is \+ small in these units:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "(G *M_s*M_e/AU^2)/unitForce;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 72 "We c an find the unit for the gravitational constant from its dimensions:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "unitG:=unitForce*(unitLen gth)^2/unitMass^2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 229 "Now the ce ntrifugal force in SI units: we assume a point at which the velocity a nd position are at right angles, so that the magnitude of angular mome ntum is given by the product of magnitudes of position and velocity ti mes mass:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "L0:=(M_e*AU*v_ e);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "L0/(unitMass*unitVel *unitLength); # angular momentum in our units:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 37 "subs(r=AU,-diff(L0^2/(2*M_e*r^2),r));" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 41 "The gravitational attraction in SI units:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "-G*M_s*M_e/AU^2; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 380 "The balance between gravitat ional attraction and the centrifugal force is not perfect, it is at th e level of the number of significant digits carried, which measn that \+ it implies that Earth's orbit is not a perfect circle, and the AU aver ages over an ellipse with eccentricity close to 0. Earth's eccentricit y is listed as 0.017. Mars has an eccentricity of 0.092, Mercury has 0 .2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 162 "No w let us figure out the balancing in the unit system introduced. We ca lculate the centrifugal force from differentiating the centrifugal pot ential expressed as " }}{PARA 0 "" 0 "" {TEXT -1 14 "L^2/(2*m*r^2):" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "subs(r=AU,diff(L0^2/(2*M_e *r^2),r))*(unitMass*unitLength^3/(unitMass*unitLength*unitVel)^2);" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "(-G*M_s*M_e/AU^2)*(unitLeng th^2/(unitMass^2*unitG));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "G/unitG;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 91 "Now, this establishe s what value one should use for G*M_s when working in this unit system :" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 148 "We \+ take a fresh start, and observe how the Earth (an object of mass m) mo ves about a fixed gravitational center using a Cartesian coordinate sy stem." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "restart; GM:=39.17 : m:=6*10^(24)/(2*10^(30)):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 66 "Th e inverse power for the potential (here we can play if we wish):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "n:=10/10;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "V:=-GM*m/r^n;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 186 "From the conservation of the orientation of the angular \+ momentum vector we know that the motion is confined to a plane, i.e., \+ the problem is two-dimensional. We label the coordinates as " }{TEXT 257 1 "x" }{TEXT -1 5 " and " }{TEXT 256 1 "y" }{TEXT -1 1 "." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "V:=subs(r=sqrt(x^2+y^2),V); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 35 "The force is calculated as fo llows:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Fx:=-diff(V,x);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Fy:=-diff(V,y);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "Note that both components of the f orce depend on both " }{TEXT 259 1 "x" }{TEXT -1 5 " and " }{TEXT 258 1 "y" }{TEXT -1 77 ". We make this more obvious by unapplying the two \+ expressions into functions:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "FX:=unapply(Fx,x,y): FY:=unapply(Fy,x,y):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 158 "We can now evaluate the force vector at some location . Let us demonstrate that the force is central, i.e., points towards t he origin of the coordinate system:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "F0:=[FX(1,2),FY(1,2)];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 63 "The y-component of the force vector at this location give n by " }{TEXT 260 1 "r" }{TEXT -1 133 " = [1,2] is twice as large as \+ the x-component, and this is sufficient for the force vector to be pro portional to the position vector." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 70 "The Newton equation leads to a pair of coupled second-order equati ons:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "NE1:=m*diff(x(t),t$ 2)=FX(x(t),y(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "NE2:=m *diff(y(t),t$2)=FY(x(t),y(t));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "Let us define some intial condition: we pick a point on the positive \+ " }{TEXT 261 1 "x" }{TEXT -1 79 "-axis, with some velocity in the y-di rection, but no velocity component in the " }{TEXT 262 1 "x" }{TEXT -1 80 "-direction, i.e., the initial position and velocity vectors are at right angles:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "x0:=1: y0:=0: vx0:=0: vy0:=6.28: # these conditions illustrate Earth's orbit " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "x0:=1: y0:=0: vx0:=0: v y0:=4.: # these conditions give substantial ellipticity" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "IC:=x(0)=x0,D(x)(0)=vx0,y(0)=y0,D(y )(0)=vy0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "sol:=dsolve(\{ NE1,NE2,IC\},\{x(t),y(t)\},numeric,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs(sol,y(t) ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "[X(1),Y(1)];" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 95 "This confirms that a trajectory is being calculated. So let us produce a graph. We need to put " }{TEXT 19 4 "X(t)" }{TEXT -1 5 " and " }{TEXT 19 4 "Y(t)" }{TEXT -1 74 " in q uotes, because they evaluate to numbers only for numerical values of \+ " }{TEXT 19 1 "t" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "plot(['X(t)','Y(t)'],t=0..3, title=\"x(t) and y(t)\",color=[blue,green],thickness=3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 36 "OK, now we plot this parametrically:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 95 "plot(['X(t)','Y(t)',t=0..3], title=\"parametric plot\",color=red,thickness=3,scaling=constrained); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT 263 11 "Exercise 1:" }}{PARA 0 "" 0 "" {TEXT -1 176 "Change the intial conditions for the equations and pr oduce the graphs. Watch how the period of the motion changes with the \+ four parameters that determine the initial condition." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 264 11 "Exercise 2:" }}{PARA 0 "" 0 "" {TEXT -1 81 "Change the potential function to a slightly dif ferent distance dependence from 1/" }{TEXT 265 1 "r" }{TEXT -1 37 ". O bserve what happens to the motion." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 98 "Let us visualize better what is g oing on by showing the vectorial quantities along the trajectory." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "Np:=30; T:=0.5;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "Lr:=[seq([X(T/Np*i),Y(T/Np*i)],i=1. .Np)]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "PLT:=plot(Lr,style=point): d isplay(PLT,scaling=constrained);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "VX:=subs(sol,diff(x(t),t)): VY:=subs(sol,diff(y(t),t) ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "with(plottools):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "scale_v:=0.1: scale_F:=300: \+ # trial and error are required to adjust these" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 36 "for i from 1 to Np do: t_i:=i*T/Np: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 105 "PLv[i] := arrow([X(t_i),Y(t_i)], evalm(s cale_v*vector([VX(t_i),VY(t_i)])), .05, .2, .1, color=green):\nod:" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "for i from 1 to Np do: t_i: =i*T/Np: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 123 "PLF[i] := arrow([X(t_ i),Y(t_i)], evalm(scale_F*vector([FX(X(t_i),Y(t_i)),FY(X(t_i),Y(t_i))] )), .05, .2, .1, color=red):\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "for i from 1 to Np do: t_i:=i*T/Np: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "PLr[i] := arrow([0,0],vector([X(t_i),Y(t_i)]), . 05, .2, .1, color=blue):\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "display(seq(display(PLT,PLv[i],PLF[i],PLr[i]),i=1..Np),insequen ce=true,scaling=constrained,view=[-1..1,-1..1]);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 92 "Let us demonstrate what happens to the potential a nd kinetic energies as a function of time:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 32 "T_kin:=t->m/2*(VX(t)^2+VY(t)^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "T_kin(0),T_kin(0.25),T_kin(0.5);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "V_pot:=t->subs(x=X(t),y=Y(t) ,V);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "E_tot:=t->T_kin(t)+ V_pot(t);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "E_tot(0),E_tot (0.25),E_tot(0.5);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "plot( ['T_kin(t)','V_pot(t)','E_tot(t)'],t=0..1,color=[red,blue,green],thick ness=3,axes=boxed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 160 "We observ e that the total energy is conserved, and that the kinetic and potenti al energies as a function of time have to be related in some way to ac hieve this." }}{PARA 0 "" 0 "" {TEXT -1 163 "When the planet is close \+ to the center of attraction it moves faster (has higher kinetic energy ) on account of the fact that the potential energy has been lowered." }}{PARA 0 "" 0 "" {TEXT -1 218 "One can establish a relationship betwe en the average values over a period (this is called the virial theorem in classical mechanics, it's details depend on the force law, and it \+ has a counterpart in quantum mechanics)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 59 "We lift an approximate value for t he period from the graph:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Period:=0.5;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "T_avg:= evalf(Int(T_kin,0..Period,5,_NCrule))/Period;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "E_avg:=evalf(Int(E_tot,0..Period,5,_NCrule))/Per iod;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "V_avg:=evalf(Int(V_ pot,0..Period,5,_NCrule))/Period;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "V_avg/E_avg;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 269 11 "Exercise 3:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 167 "Carry out these calculations for different initial conditions \+ and verify that the relationship between these average values is a fea ture of Newton's law of gravitation" }}{PARA 0 "" 0 "" {TEXT -1 144 "( inverse-squared force law). Ensure a correct choice of period by looki ng at the graph. Explore what happens for slightly different power law s, " }}{PARA 0 "" 0 "" {TEXT -1 8 "such as " }{TEXT 271 1 "n" }{TEXT -1 11 "=12/10, or " }{TEXT 270 1 "n" }{TEXT -1 11 "=8/10, etc." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 257 "" 0 " " {TEXT 272 28 "2) The radial motion problem" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 23 "The radial problem for " } {TEXT 273 1 "r" }{TEXT -1 1 "(" }{TEXT 274 1 "t" }{TEXT -1 126 ") is d erived from angular momentum conservation. It allows us to determine t he trajectory from a single differential equation." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 40 "First we calculate the an gular momentum:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(lin alg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "L0:=evalm(m*crossp rod([x0,y0,0],[vx0,vy0,0]));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "t0:=2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "evalm(m*cros sprod([X(t0),Y(t0),0],[VX(t0),VY(t0),0]));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 235 "We observe that the calculated trajectory does observe t he conservation of angular momentum approximately. This is a test if t he numerical accuracy of the algorithm. We also see that the angular m omentum vector is perpendicular to the " }{TEXT 267 1 "x" }{TEXT -1 1 "-" }{TEXT 266 1 "y" }{TEXT -1 7 " plane." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 112 "The Euclidean norm is called the \+ Frobenius norm in the linalg package (square-root of sum of elements s quared). " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "L_mag:=norm(L0 ,frobenius);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "The one-dimensio nal problem of solving for the distance from the center of attraction \+ as a function of time:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "V r:=-GM*m/r^n;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "Veff:=Vr+L _mag^2/(2*m*r^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "Feff:= unapply(-diff(Veff,r),r);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "DE:=m*diff(r(t),t$2)=Feff(r(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "r0:=sqrt(x0^2+y0^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "vr0:=0; # this is for our original choice of perpendi cular position and velocity vectors." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "IC_r:=r(0)=r0,D(r)(0)=vr0:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "#sol_r:=dsolve(\{DE,IC_r\},r(t)); # doesn't work w ell" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "sol_r:=dsolve(DE,r(t ));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "sol_r:=dsolve(\{DE,I C_r\},r(t),numeric,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 20 "R:=subs(sol_r,r(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "plot('R(t)',t=0..2,thickness=3);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 102 "We have almost balanced the centrifugal against t he gravitational forces in the case of Earth's orbit." }}{PARA 0 "" 0 "" {TEXT -1 83 "For the other initial condition a substantial variatio n in radial distance results." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 100 "We proceed with a discussion of t he stability of circular orbits for different power-law potentials." } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 127 "We choo se simpler units by picking an attractive potential of form -k/r^n, an d use a mass of unit 1 to move in this force field" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 101 "n:=3/2; # t ry this with n=3 and observe; then try n=2.0001; also observe what hap pens for n=1.999x..." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 125 "The stab ility condition: the centripetal acceleration needed for uniform circu lar motion is provided by the attractive force." }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 33 "STC:=0=Lsq/(m*r0^3)-n*k/r0^(n+1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Lsq:=solve(STC,Lsq);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "V_eff:=Lsq/(2*m*r^2)-k/r^n;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "F_eff:=unapply(diff(-V_eff,r ),r);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "NE:=m*diff(r(t),t$ 2)=F_eff(r(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "m:=1; k: =1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "r0:=1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "IC:=r(0)=r0,D(r)(0)=-1/100;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "sol:=dsolve(\{NE,IC\},r(t),n umeric,output=listprocedure):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "R:=subs(sol,r(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "plot('R(t)',t=0..100,thickness=3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 173 "A harmonic solution which drives the particle back towards osc illations about the circle with radius r0, when circular orbits are st able. Otherwise we get runaway solutions." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 109 "There is a simple graphical answe r to the question as to why the problem of circular orbits becomes uns table:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "P1:=plot(1/(2*r^2 )-1/r,r=0..10,color=red):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "P2:=plot(1/(2*r^2)-1/r^3,r=0..10,color=blue):" }}}{EXCHG {PARA 0 " > " 0 "" {MPLTEXT 1 0 80 "plots[display](P1,P2,view=[0..5,-2..1],thick ness=3,title=\"Effective Potential\");" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "The blue curve corresponds to the unstable case. We are \+ also perhaps interested in a graph of the radial force." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "P3:=plot(-diff(1/(2*r^2)-1/r,r),r=0 ..10,color=red,numpoints=300):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "P4:=plot(-diff(1/(2*r^2)-1/r^3,r),r=0..10,color=blue,numpoints =300):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "plots[display](P3 ,P4,view=[0..5,-1..0.5],thickness=3,title=\"Effective Force\");" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "The red curve represents a restori ng force with equilibrium at r=1." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 150 "Another situation that works, i.e., whic h has stable circular orbits is the case of confining (linear or highe r-power) potentials in three dimensions:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "P1:=plot(1/(2*r^2)+r,r=0..10,color=red):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "P2:=plot(1/(2*r^2)+r^2,r=0..10,colo r=blue):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "plots[display]( P1,P2,view=[0..5,0..10],thickness=3,title=\"Effective Potential\");" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "P3:=plot(-diff(1/(2*r^2)+r ,r),r=0..10,color=red,numpoints=300):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "P4:=plot(-diff(1/(2*r^2)+r^2,r),r=0..10,color=blue,nu mpoints=300):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "plots[disp lay](P3,P4,view=[0..5,-10..10],thickness=3,title=\"Effective Force\"); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT 268 25 "3) The scattering problem" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 257 "Let us \+ return to the Kepler problem, and consider positive-energy solutions. \+ These correspond to orbits that pass by the center of attraction once \+ only (true comets). These solutions are also relevant for charged part icle scattering (Rutherford scattering)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 141 "Suppose a particle is very far aw ay from the scattering center so that we can ignore the attraction. Th e 'free' particle is characterized by:" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 62 "the asymptotic velocity vector, and a ssociated kinetic energy;" }}{PARA 0 "" 0 "" {TEXT -1 118 "position ve ctor, which together with the velocity vector and the particle mass de termines the angular momentum vector." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 73 "To illustrate the idea we use Maple t o generate some numerical solutions." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "restart; with(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "m:=1: r0:=[-100,5,0]: v0:=[1,0,0]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "L0:=evalm(m*crossprod(r0,v0));" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Lsq:=dotprod(L0,L0);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "n:=1;" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 10 "V:=-1/r^n;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "V:=subs(r=sqrt(x^2+y^2),V);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "T0:=m/2*dotprod(v0,v0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "E0:=evalf(T0+subs(x=r0[1],y=r0[2],V));" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 249 "We observe that the particle is n ot entirely 'free' at a distance of 100 units. In fact, the inverse-di stance-squared force law does not fall quickly enough at large distanc es to really allow the particle to be ever truly free at finite separa tions." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "F_x:=unapply(diff (-V,x),x,y);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "F_y:=unappl y(diff(-V,y),x,y):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "NE1:= m*diff(x(t),t$2)=F_x(x(t),y(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "NE2:=m*diff(y(t),t$2)=F_y(x(t),y(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "IC:=x(0)=r0[1],D(x)(0)=v0[1],y(0)=r 0[2],D(y)(0)=v0[2];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "sol: =dsolve(\{NE1,NE2,IC\},\{x(t),y(t)\},numeric,output=listprocedure):" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "X:=subs(sol,x(t)): Y:=subs (sol,y(t)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "PL0:=plot([[ 0,0]],style=point,symbol=cross,color=red,symbolsize=20):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "PL1:=plot(['X(t)','Y(t)',t=0..200], axes=boxed,color=green,thickness=3): plots[display](PL1,PL0);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 128 "PL1:=plot([seq([X(4*i),Y(4* i)],i=1..50)],style=point,symbol=diamond,axes=boxed,color=blue,thickne ss=3): plots[display](PL1,PL0);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "VX:=subs(sol,diff(x(t),t)): VY:=subs(sol,diff(y(t),t)):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "L0,evalm(m*crossprod([X(200) ,Y(200),0],[VX(200),VY(200),0]));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "Again, angular momentum is conserved." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 140 "Don't let the graph fool you: ene rgy is conserved. The spacing between the points increases, because th e scaling on the y-axis is different." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "T_kin:=t->m/2*(VX(t)^2+ VY(t)^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "T_kin(0),T_kin (200);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "V_pot:=t->subs(x= X(t),y=Y(t),V);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "E_tot:=t ->T_kin(t)+V_pot(t);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "E_t ot(0),E_tot(100),E_tot(200);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "plot(['T_kin(t)','V_pot(t)','E_tot(t)'],t=0..200,color=[red,blue ,green],thickness=3,axes=boxed);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}} {MARK "169 0 0" 103 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }