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 $2\pi$. 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 ($0\leq \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
\qquad\rightarrow \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.5\alpha)$ 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}(-2\dot P N)^{1} -\frac{1}{8}(-2 \dot P N)^{2}
+\frac{1}{16}(-2\dot P N)^{3}+...
$$
$$
[1+ (-2\dot 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 t\approx 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 t\approx 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 ( C\dot{P}t + CP_0 )\right]_0^{T1}
$$
hence
$$
\dot{P} =
ln ( C\dot{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}$.
