The interaction of magnetohydrodynamic (MHD) waves in a non-uniform, time-dependent background plasma flow is investigated using Lagrangian field theory methods. The analysis uses Lagrangian maps, in which the position of the fluid element ${\bf x}^*$ is expressed as a vector sum of the position vector x of the background plasma fluid element plus a Lagrangian displacement $\xi({\bf x},t)$ due to the waves. Linear, non-Wentzel–Kramer–Brillouin (WKB) wave interaction equations are obtained by expansion of the Lagrangian out to second order in $\xi$ and $\Delta$S, where $\Delta$S is the Lagrangian entropy perturbation. The characteristic manifolds of the waves are determined by consideration of the Cauchy problem for the wave interaction equations. The manifolds correspond to the usual MHD waves modes, namely the Alfvén waves, the fast and slow magnetoacoustic waves and the entropy wave. The relationships between the characteristic manifolds, and the ray equations of geometrical MHD optics are developed using the theory of Cauchy characteristics for first-order partial differential equations. The first-order differential equations describing the singular manifolds are the dispersion equations for the MHD eigenmodes, where the wave vector ${\bf k}=\nabla\phi$ and frequency $\omega=-\phi_t$ correspond to the characteristic manifolds $\phi({\bf x},t)=$ constant. The form of the characteristic manifolds for both time-dependent and steady MHD flows are developed. The bi-characteristics for steady MHD waves in a steady background flow are related to the group velocity surface and Mach cone for the waves, and determine when the flow is elliptic, hyperbolic, or of mixed hyperbolic-elliptic type. The wave interaction equations are decomposed into coupled equations for the compressible and incompressible perturbations.