We occasionally get asked questions about methods of Perturbation Theory in Maple, including the Lindstedt-Poincaré Method. Presented here is the most famous application of this method.

**Introduction**

During the dawn of the 20th Century, one problem that bothered astronomers and astrophysicists was the precession of the perihelion of Mercury. Even when considering the gravity from other planets and objects in the solar system, the equations from Newtonian Mechanics could not account for the discrepancy between the observed and predicted precession.

One of the early successes of Einstein's General Theory of Relativity was that the new model was able to capture the precession of Mercury, in addition to the orbits of all the other planets. The Einsteinian model, when applied to the orbit of Mercury, was in fact a non-negligible perturbation of the old model. In this post, we show how to use Maple to compute the perturbation, and derive the formula for calculating the precession.

In polar coordinates, the Einsteinian model can be written in the following form, where `u(theta)=a(1-e^2)/r(theta)`

, with `theta`

being the polar angle, `r(theta)`

being the radial distance, `a`

being the semi-major axis length, and `e`

being the eccentricity of the orbit:

```
# Original system.
f := (u,epsilon) -> -1 - epsilon * u^2;
omega := 1;
u0, du0 := 1 + e, 0;
de1 := diff( u(theta), theta, theta ) + omega^2 * u(theta) + f( u(theta), epsilon );
ic1 := u(0) = u0, D(u)(0) = du0;
```

The small parameter `epsilon`

(along with the amount of precession) can be found in terms of the physical constants, but for now we leave it arbitrary:

# Parameters.
P := [
a = 5.7909050e10 * Unit(m),
c = 2.99792458e8 * Unit(m/s),
e = 0.205630,
G = 6.6740831e-11 * Unit(N*m^2/kg^2),
M = 1.9885e30 * Unit(kg),
alpha = 206264.8062,
beta = 415.2030758
];
epsilon = simplify( eval( 3 * G * M / a / ( 1 - e^2 ) / c^2, P ) ); # approximately 7.987552635e-8

Note that `c`

is the speed of light, `G`

is the gravitational constant, `M`

is the mass of the sun, `alpha`

is the number of arcseconds per radian, and `beta`

is the number of revolutions per century for Mercury.

We will show that the radial distance, predicted by Einstein's model, is close to that for an ellipse, as predicted by Newton's model, but the perturbation accounts for the precession (42.98 arcseconds/century). During one revolution, the precession can be determined to be approximately

`sigma = simplify( eval( 6 * Pi * G * M / a / ( 1 - e^2 ) / c^2, P ) ); # approximately 5.018727337e-7`

and so, per century, it is `alpha*beta*sigma`

, which is **approximately 42.98 arcseconds/century**.

It is worth checking out this question on Stack Exchange, which includes an animation generated numerically using Maple for a similar problem that has a more pronounced precession.

**Lindstedt-Poincaré Method in Maple**

In order to obtain a perturbation solution to the perturbed differential equation `u'+omega^2*u=1+epsilon*u^2`

which is periodic, we need to write both `u`

and `omega`

as a series in the small parameter `epsilon`

. This is because otherwise, the solution would have unbounded oscillatory terms ("secular terms"). Using this Lindstedt-Poincaré Method, we substitute arbitrary series in `epsilon`

for `u`

and `omega`

into the initial value problem, and then choose the coefficient constants/functions so that both the initial value problem is satisfied and there are no secular terms. Note that a first-order approximation provides plenty of agreement with the measured precession, but higher-order approximations can be obtained.

To perform this in Maple, we can do the following:

```
# Transformed system, with the new independent variable being the original times a series in epsilon.
de2 := op( PDEtools:-dchange( { theta = phi/b }, { de1 }, { phi }, params = { b, epsilon }, simplify = true ) );
ic2 := ic1;
```

`# Order and `