METHOD 3 Predicting Periastron Times and Cycle Durations from Orbital Phase ($Phi$)
In the answer by Stan Liou he uses a Taylor Series approximation of the Mean Anomaly to derive a nice formula which determines the CPTS (Cumulative Perihelion Time Shift) value as a function of $t^2$. This formula produces results very close to those graphed by Weisberg & Taylor. As it took me a while to understand how the Mean Anomaly can be applied for a decaying-period orbit I thought it useful to note here what I learnt and present a slightly different way of obtaining a formula to predict periastron times.
Mean Anomaly basically indicates the phase of the orbit at a particular epoch, e.g. what time fraction of the orbit period has been completed. For Mean Anomaly the fraction is scaled by $2pi$. Let $Phi(t)$ be the phase of the orbit at time $t$ such that at the start of the orbit $t=0$ and $Phi(0)$ = 0 and at the end of the orbit when $t=P$ (the period of orbit) $Phi(P)$ = 1. I will assume that an orbit starts at one periapsis and finishes at the next periapsis.
In a decaying-period orbital system the value of $P(t)$ changes with $t$. Here $P$ is the period of orbit. In a decaying-period orbit I will assume that orbital phase $Phi(t)$ is completely specified by time $t$ and either $P(t)$ or $F(t)$, i.e. aspidial precession, progressive changes in semi-major axis length or eccentricity do not lead to significant additional changes in phase.
In a decaying-period system the concept of period at a given epoch is somewhat abstract, it might be construed as a hypothetical period that would occur thereafter if the force causing period change ceased at that particular moment. Instead of period we can think in terms of $F(t)=1/P(t)$ where $F$ is the orbital frequency.
I like to think of $F(t)$ as "the rate of change of phase with time", i.e. $F(t) = dot{Phi}(t)$. Picture an imaginary phase clock comprising a circular dial whose circumference has markings running from $Phi=0$ around in a full circle to $Phi=1$. The phase $Phi(t)$ of our subject system at any given epoch $t$ will be indicated by a marker at a particular point on the circumference of the dial. Then $F(t)$ can be thought of as the speed at which the marker is moving around the circumference of the dial at a given epoch $t$. In a particular short interval of time $delta t$ the (system state) marker will move a certain distance around the dial and this distance travelled will indicate the change in phase. The "travel" (change of phase) will depend on the value of $F$ during that time as per $delta Phi approx delta t * F(t)$. This is only approximate because F(t) changes during the interval $delta t$.
Now let us assume that the time interval $T$ between succesive periaspes is known and that the time of the first periapsis is at $t=0$. If $P(t)$ (and therefore $F(t)$ too) changed in a step-wise manner between time intervals, then we could write
$$
sum_{i=1}^{i=N=T/delta t} frac{delta t}{P(t)}
=
sum_{i=1}^{i=N=T/delta t} delta t F(t)
=
sum_{i=1}^{i=N=T/delta t} delta Phi(t)
= 1.
$$
To represent continuously-changing $P$ and $F$ we reduce $Delta t$ towards zero and obtain these integral functions
$$
int_{t=0}^{t=T} frac{1}{P(t)}, mathrm{d}t
=
int_{t=0}^{t=T} F(t), mathrm{d}t
=
int_{t=0}^{t=T} dot{Phi}(t), mathrm{d}t
= 1.
$$
At any time $tau$ during the orbit ($0leq tau leq 1$) the current instantaneous phase is given by
$$
Phi(tau) =
int_{t=0}^{t=tau} dot{Phi}(t), mathrm{d}t
=
int_{t=0}^{t=tau} F(t), mathrm{d}t
=
int_{t=0}^{t=tau} frac{1}{P(t)}, mathrm{d}t.
$$
and given a constant value of $dot P$ (the rate of change of period), and $P_o=P(0)$ we obtain
$$
Phi(tau)
=
int_{t=0}^{t=tau} frac{1}{P_o + dot P t}, mathrm{d}t
qquad mathrm{and} qquad
dot{Phi}(tau) = F(tau) = frac{1}{P(tau)} = frac{1}{P_o + dot P tau}.
$$
Using
$$
frac{mathrm{d}}{mathrm{d}t} left(
frac{1}{f(x)}
right)
=
frac{-dot f(x)}{f(x)^2}
qquad mathrm{and} qquad
f(x) = P_o + tau dot P
qquad mathrm{and} qquad
dot f(x) = dot P
$$
we get
$$
ddot Phi(tau)
=
frac{mathrm{d}}{mathrm{d}t} left(
frac{1}{P_o + tau dot P}
right)
=
frac{-dot P}{P_o^2 + (tau dot P)^2 + 2 P_o tau dot P}
$$
and for $tau = 0$ we get
$$
Phi(0) = 0
qquad mathrm{and} qquad
dot{Phi}(0) = frac{1}{P_o }
qquad mathrm{and} qquad
ddot Phi(0)
= frac{-dot P}{P_o^2}.
$$
PREDICTING PERIASTRON TIMES
We can obtain an expression for CPTS by copying Stan Liou's approach , but using phase ($Phi$) instead of Mean Anomaly ($M$).
Write the phase at time $tau$ as a Maclaurin series (a Taylor series with $a=0$):
$$begin{eqnarray*}
Phi(tau) equiv int_{t=0}^{t=tau} dot{Phi}(t),mathrm{d}t
&=& Phi(0) + dot{Phi}(0) tau + frac{1}{2}ddot{Phi}(0) tau^2 + mathcal{O}(tau^3)\
&=& 0 +frac{1}{P_0}tau - frac{dot{P}}{2P_o^2}tau^2+mathcal{O}(tau^3)text{.}
end{eqnarray*}$$
Dropping the third order and residual error terms (Stan Liou's answer explains justification for this) we get
$$
Phi(tau) approx frac{1}{P_0}tau- frac{dot{P}}{2P_o^2}tau^2.
$$
Let $N$ count the number of completed orbits at some time $tau_N$. Then the phase at the end of orbit number $N$ is given by $Phi(tau_N) = N$. In a steady system (A) the orbital period $P_o$ stays constant $dot P=0$ and so the time $tau_{NA}$ at which the $N$th orbit completed is given exactly by $tau_{NA} = NP_o$ and thus $N=tau_{NA}/P_o$. In a decaying-period system (B) the Nth orbit completes at time $tau_{NB}$ when $Phi(tau_{NB}) = N$ but with $Phi(tau_{NB})$ depending on a non-zero value of $dot P$.
$$
begin{align}
&Phi(tau_{NA}) = N = frac{1}{P_0}tau_{NA}\
&Phi(tau_{NB}) = N approx frac{1}{P_0}tau_{NB}- frac{dot{P}}{2P_o^2}tau_{NB}^2.
end{align}
$$
We can obtain a formula for $Delta t_{N}$ which is the difference in time between $tau_{NB}$ (the epoch of the decaying system $N$th periastron) and $tau_{NA}$ (the epoch of the steady system $N$th periastron) as follows. Note that following Weisberg & Taylor we define $Delta t_{N},(=$ CPTS $)=tau_{NB} - tau_{NA}$. The $N$th periastron in the decaying period system occurs earlier than the $N$th periastron in the steady system. Therefore $Delta t_{N}$ (=$ CPTS $) will become increasingly negative as time passes (as $t$ increases).
$$
frac{tau_{NA}}{P_o} = N approx frac{tau_{NB}}{P_0}- frac{dot P}{2P_o^2} ,{tau_{NB}}^2
$$
$$
{tau_{NA}} approx {tau_{NB}} - frac{dot P}{2P_o} ,{tau_{NB}}^2
$$
$$
Delta t_N
= tau_{NB} - tau_{NA} approx + frac{dot P}{2P_o} ,{tau_{NB}}^2
$$
which is Stan Liou's approximation for $Delta t_N$.
So what value is this equation?
In general assume that we have determined the time $t_0$ of the $0$th periastron and measured accurately an initial orbital period $P_o$ and we keep track of the observed periastra.
Case 1 - we observe that the $N$th periastron occurs at a particular time ($tau_{NB}$). Using $tau_{NA}=N P_o$, we can easily calculate $Delta t_N = tau_{NB} - tau_{NA}$. Then, using Stan Liou's approximation of $Delta t_N$ we can obtain an empirical estimate of $dot P$ from
$$
{dot P} = frac{2 P_o ,Delta t_N} {tau_{NB}}^2.
$$
Case 2 - we have a theoretical formula which, given our measurement of $P_o$ predicts the value of $dot P$. Given various values of $N$ and using $tau_{NA}=N P_o$ we can easily calculate the hypothetical epochs ($tau_{NA}$ for various $N$) of the hypothetical steady-system periastra. Now we wish to plot a curve showing the theoretical values of $Delta t_N$ at various hypothetical periastron epochs ($tau_{NA}$). But Stan Liou's equation for $Delta t_N$ uses $tau_{NB}$ not $tau_{NA}$. We need to find $Delta t_N$ as a function of $tau_{NA}$. Therefore we would need to solve the following quadratic equation for $Delta t_N$ in terms of $tau_{NA}$ and $K$, where $K=frac{dot P}{2 P_o}$:
$$
Delta t_{N} = K {tau_{NB}}^2 = K ({tau_{NA}}^2 + 2 tau_{NA}Delta t_{N} +{Delta t_{N}}^2)
$$
$$
0 = (K{tau_{NA}}^2 + 2K tau_{NA}Delta t_{N} +K{Delta t_{N}}^2) -Delta t_{N}
$$
$$
0 = (K){Delta t_{N}}^2 +(2K tau_{NA} -1)Delta t_{N} +(K{tau_{NA}}^2)
.$$
The following way of getting an equation for $Delta t$ is more protracted but it will provide a formula for $Delta t$ as a function of $N$ (which can easily be converted into a function of $tau_{NA}$).
Starting with
$$
frac{tau_{NA}}{P_o} =frac{1}{P_0}tau_{NB}- frac{dot P}{2P_o^2} ,{tau_{NB}}^2
$$
$$
tau_{NA} =tau_{NB}- frac{dot P}{2P_o}{tau_{NB}}^2
qquadrightarrow qquad
0 = -tau_{NA} + tau_{NB} - frac{dot P}{2P_o} ,{tau_{NB}}^2
$$
Using $ x = frac{-b pm sqrt{b^2-4ac}}{2a}$ with $a=- frac{dot P}{2P_o} , b=1, c=-tau_{NA}$ we obtain
$$
tau_{NB} = frac{-(1) pm sqrt{(1)^2-4(- frac{dot P}{2P_o}).(-tau_{NA})}}{2(- frac{dot P}{2P_o})}
$$
$$
tau_{NB} = frac
{1 pm sqrt {1 - 2 frac{ dot P }{P_o}.tau_{NA} }}
{dot P /{P_o}}
qquad longrightarrow qquad
tau_{NB} =
frac{P_o}{ dot P } left( 1 pm sqrt {1 - 2 frac{dot P }{P_o}.tau_{NA} },, right)
$$
Now $tau_{NA} = P_o ,N$ so we can write
$$
tau_{NB} =
frac{P_o}{ dot P } left(1 pm sqrt {1 - 2 dot P N } right)
qquad longrightarrow qquad
-Delta t = tau_{NA}-tau_{NB}
= P_o left( N -frac{1}{ dot P} left(1 pm sqrt {1 - 2 dot P N } right) right)
$$
= = = = = = =
We can try a simple approximation of this using $(1+alpha)^{0.5} = (1+ 0.5alpha)$ so that $sqrt {1 - 2 dot P N } approx 1- dot P N$ and by inspection the "$pm$" becomes a "$+$"
$$
-Delta t
= P_o left( N -frac{1}{ dot P } left(1 pm 1 - dot P N right) right)
= P_o ( N - N )
= 0
$$
This is a valid approximation but it goes a bit too far for our purposes! We know that $tau_{NA}-tau_{NB} neq 0$. So let us try expanding the term $[1 -2 dot P N]^{0.5}$ in a series expansion,
$$
[1 + (-2 dot P N)]^{0.5}
= 1 +frac{1}{2}(-2dot P N)^{1} -frac{1}{8}(-2 dot P N)^{2}
+frac{1}{16}(-2dot P N)^{3}+...
$$
$$
[1+ (-2dot P N)]^{0.5}
= 1 - dot P N -frac{1}{2}(dot P N)^{2}
-frac{1}{2}( dot P N)^{3}+...
$$
Let us limit ourselves to terms with $ dot P $ to the power of $2$ or less thus
$$
[1-2 dot P N]^{0.5} approx
1 - dot P N -frac{1}{2}dot P ^2 N^2 .
$$
We can insert this into the time difference equation
$$
-Delta t
approx P_o left( N -frac{1}{ dot P }
left(1 pm (1 - dot P N -frac{1}{2} dot P ^2N^2 right) right)
$$
by inspection the "$pm$" becomes a "$-$"
$$-Delta tapprox P_o left( N -frac{1}{ dot P }
left( + dot P N +frac{1}{2}dot P ^2 N^2 right) right)
qquad longrightarrow qquad
-Delta tapprox P_o left( N - left( N +frac{1}{2} dot P N^2) right) right)
$$
$$
-Delta t approx -frac{1}{2}P_o dot P , N^2
qquad longrightarrow qquad
Delta t approx frac{1}{2}P_o dot P , N^2
$$
From the Weisberg & Taylor paper we are given $P_o=27906.979587552 s$ and $dot P = -2.40242 *10^{-12} s/s$
so $frac{1}{2}P_o dot P = -0.5 * 27906.979587552 * 2.40242 *10^{-12} = -3.352216747 * 10^{-08} s$.
With $N= 10,000$ cycles we get $Delta t = -3.352217 s$.
Note that we can make the substitution ${tau_{NA}}^2 = N^2 Po^2$ into the equation for $Delta t$ to give
$$
Delta t approx frac{1}{2}P_o dot P , N^2 = frac{1}{2} frac {dot P}{P_o} ,{tau_{NA}}^2
.$$
The difference between this $Delta t$ and the value of $Delta t$ from Stan Liou's equation is
$$
frac{1}{2} frac {dot P}{P_o} left( {tau_{NA}}^2 -{tau_{NB}}^2 right)
$$
$$
=
frac{1}{2} frac {dot P}{P_o} left(
{{tau_{NA}}^2
-tau_{NA}^2 +2 tau_{NA} Delta t + Delta t^2} right)
$$
$$
=
frac{1}{2} frac {dot P}{P_o} left(
{2 tau_{NA} Delta t + Delta t^2} right)
$$
and the difference as a fraction of $Delta t$
$$
=
frac{1}{2} frac {dot P}{P_o} left(
{2 tau_{NA} + Delta t} right)
.$$
For the present case $(dot P/P_o approx 8.6 * 10^{-17})$ and so, after 34,000 periastron cycles (948,838,000 s, about 30 years) the difference between the two approximations of $Delta t$ would be only $0.8 * 10^{-8} s$ which is much smaller than the time resolution of the observations.
PREDICTING CYCLE DURATIONS
We can also derive an expression for the cycle duration $D$, as follows.
The start periaston is at $t=0=T0$ and the first following periaston is at $t=T1$ then
$$
Phi(T1)=int_{t=0}^{t=T1}frac{1}{P(t)},mathrm{d}t = 1
$$
this can be expressed as
$$
1 =int_{t=0}^{t=T1}frac{1}{ dot{P} t + P_0 },mathrm{d}t
= left[ frac{1}{dot{P}} ln ( Cdot{P}t + CP_0 )right]_0^{T1}
$$
hence
$$
dot{P} =
ln ( Cdot{P}T1 + CP_0 ) - ln ( CP_0 )
=
ln left( frac{ C ( dot{P}T1 + P_0 )} { C ( P_0 )} right)
=
ln left( frac{ dot{P}T1 + P_0 } {P_0 } right)
$$
$$
rightarrow exp(dot{P})
= frac{ dot{P}T1 + P_0 } {P_0 }
rightarrow
P_0 exp(dot{P}) -P_0
= dot{P}T1
$$
giving us
$$
rightarrow
T1 =
P_0 frac{exp(dot{P}) -1} {dot{P}}
$$
so the duration of the first orbit (from start periastron to first following periastron) is $D_1 = T1-T0 = T1$ thus
$$
D_1 =
D_0 frac{exp(dot{P}) -1} {dot{P}}
$$
and we can generalize this to
$$
D_{N} =
D_{N-1} left( frac{exp(dot{P}) -1} {dot{P}} right)
.$$
For example using the made-up value $dot P = -2.34,10^{-8}$ we obtain
$$
D_{N+1} =
D_{N} frac{exp(-2.34*10^{-5}) -1} {-2.34*10^{-8}}
= 0.999 999 988 D_{N}
$$
However, using standard arithmetical software (such as Excel) when we try to calculate using $dot P = -2.34,10^{-10}$ we get nonsense results because of truncation errors. The approximate value of $dot P$ reported for the Hulse-Taylor system is about $-2.34,10^{-12}$.
We can analyze the formula for $D_{N+1}$ using series expansions. A Taylor Series expansion of $exp(x)/x -1/x$ is given by WolframAlpha as
$$
exp(x)/x -1/x= 1 +
frac{x}{2}
+ frac{x^2}{6}
+ frac{x^3}{24}
+ frac{x^4}{120}
+ frac{x^5}{720}
+ frac{x^6}{5040}
+ frac{x^7}{40320}
+ frac{x^8}{362880}
+ frac{x^9}{3628800}
+ frac{x^{10}}{39916800}
+O(x^{11})
$$
or
$$1 + frac{x^{2-1}}{2} + frac{x^{3-1}}{3*2} + ... frac{x^{i-1}}{i!} + ...$$
No easilly computable expression or approximation is obvious at present.
Proceeding anyway, it is clear that the duration of any subsequent orbit $N$ can be computed from
$$
D_N
=
D_0 left(frac{exp^{dot P} -1} {dot P} right) ^N
=
D_0 left(frac{1}{dot P} right)^N left(exp^{dot P} -1 right) ^N
$$
It is interesting to compare this expression for cycle duration with that which was used as the basis of the coarse binomial solution (see other answer: Method 4)
$$D_1=D_0(1+dot P)^1$$
The ratio of durations Phase-based to Coarse binomial based is
$$
D_0 left( frac{exp(dot P)-1}{dot P} right)^{1}
:D_0(1+dot P)^1
$$
becoming
$$
frac{exp(dot P)-1}{dot P}
:1+dot P
qquad rightarrow
exp(dot P)-1
:dot P +dot P^2
qquad rightarrow
exp(dot P)
: 1+dot P +dot P^2
$$
applying the series expansion of $exp(dot P)$ we obtain
$$
qquad rightarrow
1+dot P + frac{dot P^2}{2} + frac{dot P^3}{6} + frac{dot P^4}{24} +,...,
: 1+dot P +dot P^2
$$
Clearly the two expressions give different values for cycle duration $D$ and the difference appears at the term in $dot P^2$ and seems to be no bigger than $frac{dot P^2}{2}$.