Multi-particle azimuthal cumulants are measured as a function of centrality and transverse momentum using 470 μb−1 of Pb+Pb collisions at sNN = 5.02 TeV with the ATLAS detector at the LHC. These cumulants provide information on the event-by-event fluctuations of harmonic flow coefficients vn and correlated fluctuations between two harmonics vn and vm. For the first time, a non-zero four-particle cumulant is observed for dipolar flow, v1. The four-particle cumulants for elliptic flow, v2, and triangular flow, v3, exhibit a strong centrality dependence and change sign in ultra-central collisions. This sign change is consistent with significant non-Gaussian fluctuations in v2 and v3. The four-particle cumulant for quadrangular flow, v4, is found to change sign in mid-central collisions. Correlations between two harmonics are studied with three- and four-particle mixed-harmonic cumulants, which indicate an anti-correlation between v2 and v3, and a positive correlation between v2 and v4. These correlations decrease in strength towards central collisions and either approach zero or change sign in ultra-central collisions. To investigate the possible flow fluctuations arising from intrinsic centrality or volume fluctuations, the results are compared between two different event classes used for centrality definitions. In peripheral and mid-central collisions where the cumulant signals are large, only small differences are observed. In ultra-central collisions, the differences are much larger and transverse momentum dependent. These results provide new information to disentangle flow fluctuations from the initial and final states, as well as new insights on the influence of centrality fluctuations. [Figure not available: see fulltext.].