The Moon's orbit around the Earth is not a simple Kepler ellipse. It's dynamically distorted by tidal forces, primarily from the Sun.
Newton showed that Kepler's 3 laws arise from the inverse square law of gravity. When the mass of the orbiting body is negligible compared to the primary, its motion follows Kepler's laws. This is known as a one-body system. Newton also showed that a two-body system, where the mass of the orbiting body is not negligible, can be reduced to an equivalent one-body system.
However, when there are three or more bodies, things become complicated, especially when the mass ratios are not negligible. Much of the development of celestial mechanics in the centuries after Newton involved finding mathematical techniques for dealing with these complexities. The Moon is important for navigation and tide prediction, so considerable time and energy was devoted to developing lunar theory.
Fortunately, the Sun is a lot more massive than the Earth-Moon system, so we get a reasonable first approximation by treating the motion of the EMB (Earth-Moon barycentre) around the Sun as a Kepler ellipse. But the orbit of the Earth and Moon about their barycentre is not a nice static ellipse. The shape and size of the orbit is distorted by the Sun, with the strength of the distortion depending on the distance and direction to the Sun.
In traditional celestial mechanics, orbits are specified as ellipses given by orbital elements.
For the Earth and Moon orbits, the reference plane is the ecliptic. The Moon's orbit has a mean inclination of ~$5.14°$, but its actual inclination varies from that by ~$0.16°$. It undergoes both kinds of orbital precession. The major axis precesses (apsidal precession) with a period of ~$3232$ days ($8.85$ years), relative to the equinox direction. The line of nodes precesses backwards with a period of ~$6798$ days (18.6 years). Both of these precessions are roughly linear, but both have noticeable non-linear wobbles. robjohn has a nice diagram of smooth nodal precession in this answer:
The largest perturbation of the Moon's orbit is the evection:
It arises from an approximately six-monthly periodic variation of
the eccentricity of the Moon's orbit and a libration of similar period
in the position of the Moon's perigee, caused by the action of the
Sun.
The next largest perturbation is the variation, which compresses the orbit in the direction towards the Sun and expands the orbit in the perpendicular direction.
The other major perturbations have smaller effects, but those effects are less symmetrical. By the mid-nineteenth century, it was realised that the major perturbations could be expressed in terms of a small number of angular terms (the Delaunay arguments). Equations involving 15-20 trigonometric terms are adequate for determining the lunar orbit to a precision of around 30 seconds (of time), and predicting the occurrence of eclipses. Predicting eclipse times to higher precision, and determining the tracks of solar eclipses requires considerably more terms. The lunar theory of E. W. Brown, developed around the turn of the twentieth century, used over 1400 terms. As you can imagine, computing accurate lunar positions for navigation involved considerable labour in the days before electronic computers.
These perturbations combine together in complex ways. However, at any point in time we can take a "snapshot" of the Moon's position and velocity vectors, and find an ideal Kepler orbit which matches it at that moment. This is known as an osculating orbit:
the osculating orbit of an object in space at a given moment in time
is the gravitational Kepler orbit (i.e. an elliptic or other conic
one) that it would have around its central body if perturbations were
absent. That is, it is the orbit that coincides with the
current orbital state vectors (position and velocity).
The JPL Development Ephemeris doesn't use orbital elements to do its calculations. It performs numerical integration of the equations of motion of the major masses in the Solar System, including individual masses of 343 asteroids. However, it can compute osculating elements for the bodies in its database. We can get a feel for the true orbit by seeing how the osculating elements vary over time.
The Moon's true trajectory quickly deviates from the trajectory calculated from its osculating elements. It's quite close for a day or so, but after 6 days the error in the ecliptic longitude is ~$0.5°$.
Before we examine these osculating elements, I'll digress slightly to explain how the various elements affect the orbital period. The synodic period of the Moon is defined in terms of ecliptic longitude. Specifically, the lunar phase is proportional to the Moon's elongation, the difference between the Moon's and the Sun's ecliptic longitudes. Thus the New Moon occurs when the elongation is $0°$, First Quarter is when the Moon is $90°$ ahead of the Sun, Full Moon is when the Moon is $180°$ ahead of the Sun, etc.
The Sun's apparent motion is in the ecliptic plane, so its ecliptic longitude equals its true anomaly plus the longitude of perihelion. But the Moon's orbit is inclined to the ecliptic, which modifies the speed of its ecliptic longitude slightly, making it slower at the nodes.
I should mention that the ecliptic frame precesses backwards with a period of ~26,000 years, so cycles measured relative to the equinox differ slightly from those measured relative to the stars.
The fundamental equation of an ellipse is:
$$r = \frac{p}{1 + e\cos\nu}$$
where $r$ is the orbital radius, $p=a(1-e^2)$, $\nu$ is the true anomaly, $a$ is the semi-major axis, and $e$ is the eccentricity. The semi-minor axis is $b=a\sqrt{(1-e^2)}$.
The orbital period and speed are given by Kepler's 3rd law:
$$n^2 = \left(\frac{2\pi}{P}\right)^2=\frac{\mu}{a^3}$$
where $P$ is the (anomalistic) period, the time between successive periapsis passages, $n$ is the mean motion, and $\mu$ is the gravitational parameter. For a two-body system, $\mu=G(m_1+m_2)$, where $G$ is the universal gravitational constant and $m_i$ are the masses. (As Wikipedia explains, it's very hard to measure $G$, but we know $\mu$ for major Solar System bodies to quite high precision).
Another important parameter is $\mathbf{h}$, the specific angular momentum. This vector is perpendicular to the orbit, and it's a constant of motion in a one- or two-body system. From Kepler's 2nd law, its magnitude $h$ is half the areal velocity. Thus $h=nab$. Also, $h^2=\mu p$, and $h=r^2\dot\nu$, where $\dot\nu$ is the rate of change of the true anomaly, i.e., it's the body's angular speed in the orbital plane. So if $h$ is constant, the angular speed of an orbiting body is inversely proportional to the square of the distance from the body to its primary. The Moon's specific angular momentum isn't constant, but it only varies by ~1%.
Here are some plots of the main osculating elements that relate to the angular speed. These plots cover the period between two New Moons, 2022-Jan-02 UTC and 2023-Jan-22 UTC, a span of ~13 synodic months, and ~14 of the other kinds of lunar month.
Eccentricity
Semi-major axis
Inclination
Mean motion
Sidereal orbit period
Specific angular momentum
Here's the Python / Sage script that generated those plots (except the angular momentum plot), using data from JPL Horizons. It can generate plots for a few other osculating elements, for any body that Horizons knows.
As you can see, all of these elements vary quite quickly, and the amount of variation is not constant. However, the total amount of variation is mostly fairly small. The semi-major axis varies by ~1%. The inclination varies by ~3%, but that's not so bad because the inclination is fairly small. The worst "offender", by far, is the eccentricity, which is why it's commonly blamed for the variation in the synodic month length.
The size of the Moon's mean eccentricity does make a significant contribution to the synodic (and sidereal) month length variation. However, because the Moon's orbit is precessing, changing eccentricity, and changing size, the pattern of synodic month length variation is even more complex.
The maximum variation in the synodic month length purely due to the eccentricities is around 8 hours. Here's a plot showing the month length (from Full Moon to Full Moon) covering 235 synodic months (a Metonic cycle, almost exactly 19 tropical years, and 254 tropical months), in a simplified system where the orbits have the same eccentricities as the real system, and the same tropical & synodic periods, but the orbits are coplanar, with no perturbations.
The vertical axis shows the difference in month length from the mean synodic length, in hours.
Here's the plotting script, which can also print that data in CSV format, with a few extra items: the month length (in days), and the true anomaly of the Moon at the moment of Full Moon (in degrees).
Interestingly, the variation from New Moon to New Moon in this system, is only ~±4 hours, and between First or Last Quarters it's ~±6 hours.
You may like to look at elements plots of the Earth, and of other Solar System bodies, to see how well-behaved they are, compared to the Moon.
