Finding the initial conditions of a conservation law hyperbolic system? Consider $u(\boldsymbol{x},t)$ to be a function that is known at $t=0 $ by the initial conditions
$$u(\boldsymbol{x},0)= f(\boldsymbol{x}),\quad u_t(\boldsymbol{x},0)= g(\boldsymbol{x})$$
Where $\boldsymbol{x}\in \mathbb{R}^3$ 
Let $q, r, s$ be functions such that
\begin{align*}
    &q_t= cu_x  \\
    &r_t= cu_y  \\
    &s_t= cu_z  \\
    &u_t =cq_x+ cr_y+ cs_z
\end{align*}
Then, apparently, up to a constant, we can determine $q(\boldsymbol{x},0),\ r(\boldsymbol{x},0)$ and $s(\boldsymbol{x},0)$ at the initial time, in term of $f$ and $g$ . 
But I stuck at how to find that.
The question is : Can we determine these initial values ??

\ tentative :
We can notice that the system of equations can be looked at as two simple equations if we consider $v= (q,r,s)$ :
\begin{align*}
    &v_t= c\ \nabla u \\
    &u_t =c\ \nabla \cdot v
\end{align*}
Notice also, that for the one dimensional case:
\begin{align*}
    &q_t= cu_x \\
    &u_t =cq_x
\end{align*}
We can easily derive that
\begin{align*}
    &c\ q(x,t)= \int_{x_0}^x u_t(\xi,t)d\xi \\
\text{So that }\quad  & q(x,0)= \frac{1}{c}\int_{x_0}^x u_t(\xi,0)d\xi= \frac{1}{c}\int_{x_0}^x g(\xi)d\xi
\end{align*}
And we see that indeed in the one dimensional case $q(x,0)$ is known in term of $g$ up to a constant. So now the question is to do the same for the 3 dimensional case.
 A: If you take the spatial and temporal derivatives of $q,\,r,\,s$, you'll quickly find,
\begin{align}
  q_{tx} &= cu_{xx} \\
  r_{ty} &= cu_{yy} \\
  s_{tz} &= cu_{zz}
\end{align}
And if you then take the second temporal derivative of $u$, you can see pretty easily,
$$ u_{tt}=c\left(q_{x}+r_{y}+s_{z}\right)_t=c\left(q_{tx}+r_{ty}+s_{tz}\right)=c^2\left(u_{xx}+u_{yy}+u_{zz}\right)$$
which is to say that this is really just the wave equation. In one dimension, you will find that this is solved using d'Alembert's formula,
$$ u\left(x,t\right)=\frac{1}{2}\left(f(x+ct)+f(x-ct)\right)+\frac{1}{2c}\int_{x-ct}^{x+ct}g(\xi)\,\mathrm d\xi.\tag{1}$$
In higher dimensions, this becomes known as Kirchhoff's formula (PDF link), which is a bit more complicated to derive, but more or less follows the same idea as the 1D case. You want to convert to spherical coordinate system and consider spherical means, defining $U$, $F$ and $G$ such that,
\begin{align}
  U(x,r,t)&=\frac{1}{\vert\partial B(x,r)\vert}\int_{\partial B(x,r)}u(\xi,t)\mathrm dS_\xi \\
  F(x,r,t)&=\frac{1}{\vert\partial B(x,r)\vert}\int_{\partial B(x,r)}f(\xi,t)\mathrm dS_\xi \\
  G(x,r,t)&=\frac{1}{\vert\partial B(x,r)\vert}\int_{\partial B(x,r)}g(\xi,t)\mathrm dS_\xi
\end{align}
where $B$ is the ball in $\mathbb{R}^d$ and $\partial B$ the surface. To recover $u$ from this, we simply take the limit of $r\to0$,
$$u(x,t)=\lim_{r\to0}U(x,r,t)$$
If we then let $\bar{U}=rU$ and likewise for $F$ and $G$, then $\bar{U}$ satisfies the wave equation and has a solution similar to Eq (1),
$$\bar{U}(x,r,t)=\frac{1}{2}\left[\bar{F}(t+r)+\bar{F}(t-r)\right]+\frac{1}{2}\int_{t-r}^{t+r}\bar{G}(y)\mathrm dy $$
Taking the limit gives,
$$ u(x,t)=\lim_{r\to0}\frac{\bar{U}}{r}=\frac{\partial}{\partial t}\left(\frac{t}{\vert\partial B\vert}\int_{\partial B}f(\xi)\mathrm dS_\xi\right)+\frac{t}{\vert\partial B\vert} \int_{\partial B}g(\xi)\mathrm dS_\xi.$$
With quite some more work, this can be reduced to what is called Kirchhoff's formula,
$$ u(x,t)=\frac{1}{\vert\partial B\vert}\int_{\partial B}\left[tg(\xi)+f(\xi)+\nabla f\cdot\left(\xi-x\right)\right]\,\mathrm dS_\xi$$


When solving the wave equation $\psi_{tt}=c^2\nabla^2\psi$, a common trick is to define variables $p$, $q$, $r$ and $s$ such that,
$$ p=c\psi_x,\quad q=c\psi_y,\quad r=c\psi_z,\quad s=\psi_t$$
and combine them into $u=\left[p,\,q,\,r,\,s\right]^T$ to form the evolution equation,
$$ \frac{\partial u}{\partial t}=-\nabla F\tag{A.1}$$
where $F$ is flux in terms of the variables $p,\,q,\,r,\,s$. For instance,
\begin{align}
  p_t &= c\psi_{xt}=c\psi_{tx}=cs_x \\
  q_t &= c\psi_{ty}=cs_y \\
  r_t &= c\psi_{tz}=cs_z \\
  s_t &= cp_x+cq_y+cr_z
\end{align}
implies,
$$ F_x=\left(\begin{array}{cccc} 0 & 0 & 0 & -c \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -c & 0 & 0 & 0\end{array}\right)\cdot\left(\begin{array}{c} p \\ q \\ r \\ s \end{array}\right) $$
I believe that what you are supposed to be doing is that, given the initial conditions $\psi(x,0)=f(x)$ and $\psi_t(x,0)=g(x)$, you apply these directly to the variables:
$$ p(x,0)=c\frac{\partial f}{\partial x},\;q(x,0)=c\frac{\partial f}{\partial y},\;r(x,0)=c\frac{\partial f}{\partial z},\;s(x,0)=g(x) $$
And integrate according to Eq (A.1) using Lax-Wendroff, rather than trying to apply Kirchhoff's formula to the equations.
