  We propose a method to constrain the variation of the gravitational constant
$G$ with cosmic time using gravitational-wave (GW) observations of merging
binary neutron stars. The method essentially relies on the fact that the
maximum and minimum allowed masses of neutron stars at a particular cosmic
epoch has a simple dependence on the value of $G$ at that epoch. GWs carry an
imprint of the value of $G$ at the time of the merger. Thus, if the value of
$G$ at merger is significantly different from its current value, the masses of
the neutron stars inferred from the GW observations will be inconsistent with
the theoretically allowed range. This enables us to place bounds on the
variation of $G$ between the merger epoch and the present epoch. Using the
observation of the binary neutron star system GW170817, we constrain the
fractional difference in $G$ between the merger and the current epoch to be in
the range $-1 \lesssim \Delta G/G \lesssim 8$. Assuming a monotonic variation
in $G$, this corresponds to a bound on the average rate of change of $-7 \times
10^{-9}~\mathrm{yr}^{-1} \le \dot{G}/G \le 5 \times 10^{-8}~\mathrm{yr}^{-1}$
between these epochs. Future observations will put tight constraints on the
deviation of $G$ over vast cosmological epochs not probed by other
observations.
