I was trying to plot precession of perihelion of Mercury using matlab. For this I am following a book Computational Physics by Nicholas J. Giordano and Hisao Nakanishi 2nd Edition. In that book authors used following force law predicted by general relativity:

$$ F \approx \frac{G M_m M_s}{r^2}\left(1+\frac{\alpha}{r^2}\right) $$

Where $\alpha$, they state to be expressed by speed of light, mass of sun, eccentricity of the orbit.

To plot the orbit I have used that equation and ode45, ode23 functions in matlab. My code is as follows:

% Initial Conditions
% Initial position in AU
xm0 = [0.47 0];
% Initial velocity in AU/Year
vm0 = [0 8.2];

inivec = [xm0;vm0]

% tspan for 1000 years
tspan = [0, 1000];

[t, y] = ode45(@sunmer, tspan, inivec);

% x positions of Mercury
xmer = y(:,[1:2]);
% y positions of Mercury
vmer = y(:,[3:4]);
hold on

% Position of Sun in plot
hold off
axis([-0.5, 0.5, -0.5, 0.5])

function for ode45 or ode23

function xpr = sunmer(t,x)
alpha = 0.0008;
xmer = x([1:2]);
vmer = x([3:4]);
r = norm(xmer);
rcubed = r^3;
rsqr = r^2;

% amer = -((G*Ms*x)/r^3)*(1+(alpha/r^2))
% G*Ms = 4*pi^2 in AU^3/Yr^2
amer = -(((4*(pi)^2))/rcubed)*xmer*(1+(alpha/rsqr));
xpr = [vmer;amer];

I am getting following plot when using ode45: enter image description here

If I use ode23 I get following figure: enter image description here

These are not the plots I was expecting. I can't understand what is wrong in my code. Is it my code or the ODE solvers?

  • $\begingroup$ When debugging make your problem simpler. Set $\alpha=0$. Do you obtain closed orbits in that case? $\endgroup$ Dec 5, 2019 at 0:08
  • $\begingroup$ Did you reduce the values for the absolute and relative tolerances requested in ode23s and ode45? Add opts = odeset('RelTol',1e-10,'AbsTol',1e-8); and use [t, y] = ode45(@sunmer, tspan, inivec, opts); $\endgroup$
    – GertVdE
    Jan 22, 2020 at 11:15
  • $\begingroup$ Also note that for orbit calculations, you better resort to symplectic (i.e. energy-conserving) integration methods. See for example this question and its answers $\endgroup$
    – GertVdE
    Jan 22, 2020 at 11:20


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.