PHYS 325 Discussion 6 : Perihelion Precession of Mercury
PHYS 325 Discussion 6 : Perihelion Precession of Mercury
By - Sagnik Saha, 10/13/2025.
Part (a)
Part (a)
Begin by defining the perturbative expansion of ũ, up to order λ^2.
In[]:=
ũ[ϕ_]:=ũ0[ϕ]+λ*ũ1[ϕ]+λ^2*ũ2[ϕ];
Plug into Eq.(4), and isolate the powers of λ. I am pulling everything to the LHS, to make the RHS = 0.
In[]:=
diffeq=Expand[D[ũ[ϕ],{ϕ,2}]+ũ[ϕ]-λ*ũ[ϕ]^2-1]
As seen above, the expression has terms up to O(λ^5), which we don't care about. Extract terms up to O(λ^1).
We use a Mathematica function called Coefficient.
Syntax: Take the expression 'diffeq', and return the coefficient of λ^0. This is the Binet equation you saw in lecture.
In[]:=
Coefficient[diffeq,λ,0]
In[]:=
Coefficient[diffeq,λ,1]
The expressions above are set to 0 to get Eqs.(6a) and (6b).
Part (b)
Part (b)
Define Eq.(9) for the ansatz for ũ1[ϕ].
NOTE: I am using lowercase for a,b,c. This is good practice in Mathematica, since there are some uppercase letters that are 'protected', meaning that they mean a specific thing to Mathematica and hence cannot be modified by the user.
Eg - Try ?C and ?E in a cell, and see what you get.
Eg - Try ?C and ?E in a cell, and see what you get.
In[]:=
ũ1[ϕ_]:=a+b*ϕ*Sin[ϕ]+c*Cos[2*ϕ]
In[]:=
ũ1[ϕ](*Testcall*)
Plug into the left most side of Eq.(8)
In[]:=
eq8withansatz=D[ũ1[ϕ],{ϕ,2}]+ũ1[ϕ]
Now, compare this to the rightmost side of Eq.(8). By comparing coefficients, it is easy to see that:
In[]:=
a=(1+ϵ^2/2);b=ϵ;c=-ϵ^2/6;
Now that we have a, b, c, define ũ1[ϕ] again .
In[]:=
ũ1[ϕ_]:=a+b*ϕ*Sin[ϕ]+c*Cos[2*ϕ]
In[]:=
ũ1[ϕ](*Testcall*)
Also define ũ0[ϕ] using Eq.(7)
In[]:=
ũ0[ϕ_]:=1+ϵ*Cos[ϕ]
Finally, using Eq.(5), define the full solution ũ[ϕ] to first order in λ.
In[]:=
ũ[ϕ_]:=ũ0[ϕ]+λ*ũ1[ϕ]
In[]:=
ũ[ϕ]
Part (c)
Part (c)
Take the derivative of ũ[ϕ].
In[]:=
dervu=D[ũ[ϕ],ϕ]
First, check that ϕ = 0 is a solution of ũ'[ϕ]=0
In[]:=
dervu/.ϕ->0
Then, check that ϕ = 2*Pi is NOT a solution of ũ'[ϕ]=0
In[]:=
dervu/.ϕ->2*Pi
Now, let's solve for δ. First, plug in ϕ = 2*Pi + δ
In[]:=
dervufullperiod=FullSimplify[dervu/.{ϕ->2*Pi+δ}]
Then, plug in first order Taylor expansion for sin and cos (see hint at the end of part (c).
In[]:=
dervufullperiod=dervufullperiod/.{Sin[δ]->δ,Cos[δ]->1}
Solve for δ using the expression above. Syntax: The Solve function gives something like {δ -> ...}, which we store in a variable δ using a rule.
In[]:=
δfullsol=δ/.Solve[dervufullperiod==0,δ]
Finally, we need the leading order in λ term, so we simply Taylor expand to first order in λ.
In[]:=
foshift=Series[δfullsol,{λ,0,1}]
Let's store only the first order term in a separate variable for convenience. The line below pulls out the first order term from the ‘list’ above. Just an idiosyncrasy of Mathematica.
In[]:=
foshift=Normal[foshift][[1]]
Part (d)
Part (d)
From the definitions given earlier, we know that:
λ3
2
GM
cℓ
And,
2
ℓ
2
ϵ
CAUTION : Now ' c' is the speed of light in vacuum, and not the parameter we derived earlier .
Using the expressions above, we can solve for λ in terms of a,c, and ϵ, and plug it into our expression for δ. We do this now.
First, make sure to clear any previously defined values of a,c.
In[]:=
Clear[a,c]
Now, make the replacement.
In[]:=
foshiftsimp=foshift/.{λ->3*G*M/(a*c^2*(1-ϵ^2))}
Plug in the values of G,M,c,a,ϵ. First three are constants, last two are given to us.
In[]:=
foshiftmerc=foshiftsimp/.{G->6.67*10^(-11),M->2*10^30,c->3*10^8,a->5.79*10^10,ϵ->0.2056}
This result is in units of 'radians'. Convert it to arcsec/century on paper (it’s just multiplication by a factor).
Plots
Plots
Let’s plot the unperturbed and perturbed orbits to see the precession explicitly. We plot over 8 orbits, i.e. ϕ = 0 to 16*π.
For convenience, let us define the perturbed and unperturbed orbit equations here, all in one place.
NOTE: The actual precession of Mercury is really tiny (as we already saw). As such, I will use a fictional value of λ = 0.01 (the real number is of order 10^(-8)) to see the precession clearly. I keep ϵ = 0.2056, i.e. the real eccentricity of Mercury’s orbit. Finally, I use r̃ = 1/ũ, to plot r̃ (ϕ) instead of r(ϕ), since the latter just amounts to scaling by a constant and hence, doesn’t affect the shape of the curve at all (alternatively, think of it as setting setting α = 1).
NOTE: The actual precession of Mercury is really tiny (as we already saw). As such, I will use a fictional value of λ = 0.01 (the real number is of order 10^(-8)) to see the precession clearly. I keep ϵ = 0.2056, i.e. the real eccentricity of Mercury’s orbit. Finally, I use r̃ = 1/ũ, to plot r̃ (ϕ) instead of r(ϕ), since the latter just amounts to scaling by a constant and hence, doesn’t affect the shape of the curve at all (alternatively, think of it as setting setting α = 1).
Define λ (fictional) and ϵ.
In[]:=
λ=10^(-2);ϵ=0.2056;
Unperturbed Orbit Equation
In[]:=
ũ0[ϕ_]:=1+ϵ*Cos[ϕ]r̃0[ϕ_]:=1/ũ0[ϕ]
Perturbed Orbit Equation (using Eqs.(5), (9), and (16)*)
In[]:=
ũ1[ϕ_]:=1+ϵ^2/2+ϕ*Sin[ϕ]*ϵ+Cos[2*ϕ]*(-ϵ^2/6);ũ[ϕ_]:=ũ0[ϕ]+λ*ũ1[ϕ];r̃[ϕ_]:=1/ũ[ϕ];
Now let’s make the plots. We use ParametricPlot. Since r̃ is in polar coordinates, we have to use x = r̃cos(ϕ), y = r̃sin(ϕ). Alternatively, you could just use the function PolarPlot.
In[]:=
ParametricPlot[{{r̃0[ϕ]*Cos[ϕ],r̃0[ϕ]*Sin[ϕ]}},{ϕ,0,16*Pi},PlotLabel->"Unperturbed (8 Orbits)",AxesLabel->{"x","y"}]
In[]:=
ParametricPlot[{{r̃[ϕ]*Cos[ϕ],r̃[ϕ]*Sin[ϕ]}},{ϕ,0,16*Pi},PlotStyle->AbsoluteThickness[1.0],PlotLabel->"Precession (8 Orbits)",AxesLabel->{"x","y"}]
In[]:=
ParametricPlot[{{r̃0[ϕ]*Cos[ϕ],r̃0[ϕ]*Sin[ϕ]},{r̃[ϕ]*Cos[ϕ],r̃[ϕ]*Sin[ϕ]}},{ϕ,0,16*Pi},PlotStyle->{AbsoluteThickness[1.0],AbsoluteThickness[1.0]},PlotLegends->{"Unperturbed","Perturbed"},AxesLabel->{"x","y"},PlotLabel->"Precession v/s no precession (8 Orbits)"]