I have a system of four fields $\mathbf{u}(t,x)=(\alpha,\gamma,\beta_{+},\beta_{-})$ over one spatial dimension and one time dimension. The system is expressed as an inital value problem and each field has an evolution equation which will depend on the other fields.
$$ \partial_{t}\alpha = -\frac{\alpha^{2}f(\beta_{+}+\beta_{-})}{\sqrt{\gamma}} \\ \partial_{t}\gamma = -\alpha\gamma(\beta_{+}+\beta_{-}) \\ \partial_{t}\beta_{+} = -\partial_{x}\left(\frac{\alpha\sqrt{f}\beta_{+}}{\sqrt{\gamma}}\right) \\ \partial_{t}\beta_{-} = \partial_{x}\left(\frac{\alpha\sqrt{f}\beta_{-}}{\sqrt{\gamma}}\right) $$ where, $f=\mathrm{constant}>0$. Initially, $\beta_{+}$ and $\beta_{-}$ are taken to be localized smooth wave pulses which then travel at constant speed in either direction ($\beta_{+}$ to the right and $\beta_{-}$ to the left). These perturbations result in similar pulses in $\alpha$ and $\gamma$.
For linear systems it is useful to take the solution as a Fourier mode $\mathbf{\hat{u}}_{k}\mathrm{e}^{-i\omega t}\mathrm{e}^{ikx}$ and plug it into the equations to obtain the dispersion relation $\omega(k)$. I have read somewhere that nonlinearity may also affect the amplitude of the nodes, so I am unsure of how to perform such an analysis in the nonlinear case.
From a different analysis, I know that something special happens when $f=1$. In that case, the pulses should retain their shape and travel indefinitely with no dispersion. Whenever $f\neq 1$, the solutions may lead to shocks in finite time. I have attempted to linearize the system but the value for $\omega$ I obtained after finding the eigenvalues of the 4x4 Jacobian contained an imaginary part (will lead to decay or amplification) even when $f=1$ so I don't trust it.
What is the dispersion relation for this system?