The article investigates a mathematical model of the Duffing oscillator with a variable fractional order derivative of the Riemann–Liouville type. The study of the model is carried out using a numerical scheme based on the approximation of the fractional derivative of the Riemann–Liouville type by a discrete analog—the fractional derivative of Grunwald–Letnikov. The adequacy of the numerical scheme is verified using specific examples. Using a numerical algorithm, oscillograms and phase trajectories are constructed depending on the values of the model parameters. Chaotic regimes of the Duffing fractional oscillator are investigated using the Wolf–Bennetin algorithm. The forced oscillations of the Duffing fractional oscillator are investigated using the harmonic balance method. Analytical formulas for the amplitude-frequency, phase-frequency characteristics, and also the quality factor are obtained. It is shown that the fractional Duffing oscillator possesses different modes: regular, chaotic, multi-periodic. The relationship between the order of the fractional derivative and the quality factor of the oscillatory system is established.