Weak shock solutions of system of conservation laws We consider the system of autonomous conservation laws
$$\mathbf{u_t} + (\mathbf{f(u)})_x = 0$$
with left state $\mathbf{u}_L = \mathbf{u}_0$ and right state $\mathbf{u}_R = \mathbf{u}_0 + \epsilon\delta\mathbf{u}$, where $\epsilon$ is small and $\mathbf{u}_0$ and $\delta\mathbf{u}$ are constant space vectors in $R^N$. We want to find N distinct families of weak shock solutions with right state and shock speed $$\mathbf{u}_{R, \alpha} = \mathbf{u}_0 + \epsilon\mathbf{r}_\alpha,   \qquad       s=s_\alpha$$ for $\alpha = 1, 2, …, n$ and find $\mathbf{r}_\alpha$ and $s_\alpha$.
The Rankine-Hugoniot condition for this system is $s(\mathbf{u}_L - \mathbf{u}_R) = \mathbf{f}(\mathbf{u}_L) - \mathbf{f}(\mathbf{u}_R)$. We need to linearise it in $\delta\mathbf{u}$ by taking the derivative with respect to $\epsilon$ and setting $\epsilon = 0$ to proceed, so with the given states I transformed it to $s(\epsilon\delta\mathbf{u}) = \mathbf{f}(\mathbf{u}_0) - \mathbf{f}(\mathbf{u}_0 + \epsilon\delta\mathbf{u})$.
My attempt at linearising:
$\lim_{\epsilon\to 0} \frac {d}{d\epsilon}[\mathbf{f}(\mathbf{u}_0) - \mathbf{f}(\mathbf{u}_0 + \epsilon\delta\mathbf{u})] = \mathbf{f}’(\mathbf{u}_0)-\mathbf{f}’(\mathbf{u}_0)\delta\mathbf{u}$
I feel like I’ve linearised the wrong thing as I don’t see how this helps us to find the right state and shock speed and I’m really struggling to get started. What am I doing wrong or what am I missing?
Any help would be massively appreciated!
 A: Good start. As discussed in this related thread, we can write the Rankine-Hugoniot conditions
\begin{aligned}
s_\alpha \left({\bf u}_{R,\alpha} - {\bf u}_{R,\alpha-1}\right) &= {\bf f}({\bf u}_{R,\alpha}) - {\bf f}({\bf u}_{R,\alpha-1})\\
&= {\bf f}({\bf u}_{0} + \epsilon {\bf r}_\alpha) - {\bf f}({\bf u}_{0} + \epsilon {\bf r}_{\alpha-1}) \\
= \epsilon s_\alpha \left({\bf r}_\alpha - {\bf r}_{\alpha-1}\right)
&\simeq \epsilon {\bf A}({\bf u}_{0}) \left( {\bf r}_\alpha - {\bf r}_{\alpha-1} \right)
\end{aligned}
where we have used the Taylor series approximation
$$
{\bf f}({\bf u}_{0} + \epsilon \, \delta{\bf u}) \simeq {\bf f}({\bf u}_{0}) + \epsilon {\bf A}({\bf u}_{0}) \delta{\bf u}
$$
involving the Jacobian matrix ${\bf A} = {\partial {\bf f}}/{\partial \bf u}$. Finally, $s_\alpha$ can be approximated by an eigenvalue of the Jacobian matrix ${\bf A}({\bf u}_{0})$ for small $|\epsilon| > 0$, where the difference ${\bf r}_\alpha - {\bf r}_{\alpha-1}$ belongs to the corresponding eigenspace. It remains to sort the eigenvalues of ${\bf A}({\bf u}_{0})$ by ascending order to finish (see also this related post for an illustration).
