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$ – Wolfgang Bangerth Dec 5 '19 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 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 at 11:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.