Is the shortest day duration constant?
By looking at the link provided by James K, I was able to write a Python 3 script that calculates the duration of the longest night for a predefined location (Boston, MA) for 160 years. The script uses the PyEphem module. The longest night becomes shorter by ~72 msec per year, and its length also oscillates with the period of ~19 years.



Here's the script that I used:

import ephem
import matplotlib.pyplot as plt
import matplotlib
import numpy as np

location = ephem.Observer()
location.lat, location.lon = '42.3601', '-71.0589'

years = range(1900, 2060)
longest_nights = []
for year in years:
    location.date = ephem.next_winter_solstice('{}-12-12'.format(year))
    sunset = location.next_setting(ephem.Sun())
    sunrise = location.next_rising(ephem.Sun())
    if sunrise &lt; sunset: # Look at the night before!
        location.date -= 1       
        sunset = location.next_setting(ephem.Sun())
    longest_nights.append(sunrise - sunset)

longest_nights = np.array(longest_nights)
longest_nights = (longest_nights - longest_nights.mean()) * 24 * 3600 # In seconds

# Find the slope
xp = np.poly1d(np.polyfit(years, longest_nights, 1))
np_years = np.array(years).reshape(len(years), 1)

# Plot
matplotlib.style.use("ggplot")
plt.plot(years, longest_nights, "-o")
plt.ylabel("Deviation from the mean, sec")
plt.xlabel("Year")
plt.title("Duration of the longest night in a year in Boston")
plt.plot(years, xp(years), "-")
plt.arrow(2017, plt.ylim()[0], 0, plt.ylim()[1] - plt.ylim()[0], color='k')
#plt.savefig("solstice.png")
plt.show()