Estimated solution to system of equations with phase-shifted functions Forgive my first attempt at MathJax. I have a system of $n$ equations of the form
$$
v_j(t) = \sum_{i=0}^{m-1} \frac 1 {|\vec p_i - \vec q_j|}
u_i \left(t - \frac {|\vec p_i - \vec q_j|} s \right)
$$
Variables
$ i \in \mathbb{Z}, 0 \leq i < m $
$ j \in \mathbb{Z}, 0 \leq j < n $
$ m \in \mathbb{Z}, m > 0 $, unknown, on the order of ~1-20
$ n \in \mathbb{Z}, n > 1 $, known, on the order of ~2-10
$ u_i(t) \in \Re $, unknown functions
$ v_j(t) \in \Re $, known functions
$ \vec p_i = 
   [\begin{matrix} \hat x & \hat y & \hat z \end{matrix}]
   \left[ \begin{matrix} p_{ix} \\ p_{iy} \\ p_{iz} \end{matrix} \right]
$, components $\in \Re$, unknown
$ \vec q_j = 
   [\begin{matrix} \hat x & \hat y & \hat z \end{matrix}]
   \left[ \begin{matrix} q_{jx} \\ q_{jy} \\ q_{jz} \end{matrix} \right]
$, components $\in \Re$, known
$ s \in \Re, s > 0 $, known constant
$ t \in \Re $, domain variable over which $v_j$ are known for some finite $t_0 \leq t \leq 0$
Other Constraints
Practical constraints are known for all variables. In particular, it is expected that, most likely:
$|\vec p_{i1} - \vec p_{i2}| \gtrsim |\vec q_{j1} - \vec q_{j2}|, i_1 \neq i_2, j_1 \neq j_2  $
$u_i(t) \gtrsim v_j(t)$ for any $i, j$ due to the term $\frac 1 {|\vec p_i - \vec q_j|}$
Discussion
I need to solve for all unknowns. Since there are many more degrees of freedom than there are equations, I realize that the unknowns will likely have to be found with some numerical iterative approximation method in a computer program of my design. $m$ should be minimized while still yielding an accurate estimation of all $u_i$.
One back-of-the-envelope attempt I made at a solution for $m=2, n=2$ ended up yielding a recurrence relation of the form
$$ u_0(t) + \mu_0 u_0(t - \phi_0) = v_0(t - \phi_1) + \mu_1 v_1(t - \phi_2) $$
where I assumed that this was inside of an iteration that calculated all $\mu$ and $\phi$ based on estimates for $\vec p$.
Recurrence Solution
I worked out the following generic solution to a recurrence equation of the type I (maybe) need.
Given some original $f(t)$ and $g(t) = f(t) + \mu f(t-\phi)$, the original can be restored via
$$f(t) = \sum_{i=0}^\infty (-\mu)^i g(t-\phi i)$$
This suggests that the solution to the above "envelope" equation is
$$ u_0(t) = \sum_{i=0}^\infty (-\mu_0)^i
\left(
v_0(t - \phi_1 - i \phi_0) + \mu_1 v_1(t - \phi_2 - i \phi_0)
\right) $$
Questions


*

*Did I correctly solve the recurrence relation above?

*What are some approaches I can take to solve the system of equations numerically?

*Will all cases yield recurrence relations?

*Does this problem belong to a known class of problems?

*Does it make a big difference if this problem is framed with $t$ discrete, and $u_i$ and $v_j$ discrete series instead of continuous functions? How would the solution approach change?

 A: The equation can be expressed in the general form of a convolution.
$v = G \ast u$,
where $u(x, t)$ is a distribution of sources which generate the waves (stationary deltas centred at each point source, in your case),
$G(x, t) = \frac 1 {2 \tau |x|} \delta(t - \frac {|x|} s)$
is the convolution kernel, and $v(x, t)$ is the resulting sound wave. G is the sound wave that results when a single perfect bang occurs at the origin at $t = 0$ and extends outward in a sphere. The convolution is done over all space and time. The goal now is to reverse the effect of the convolution, which can be done using the space-time Fourier transform:
$\tilde f(k, \omega) = \mathrm F(f)(k, \omega) = \tau^{-2} \int_{xt} e^{j(\omega t - k \cdot x)} f(x, t)$
Doing the Fourier transform to v gives you
$\tilde v = \tau^2 \tilde G \tilde u$.
Dividing by $\tau^2 \tilde G$,
$\tilde u = \frac {\tilde v} {\tau^2 \tilde G}$.
Then it's only a matter of finding $\tilde G$ mathematically, which shouldn't be too hard, and getting enough samples of $v$ to make a good enough approximation of $\tilde v$. I'm actually working on $\tilde G$ right now. I'm almost there. It'll have to wait until morning, though.
Since I left $u$ as a general source distribution, it could be anything, and not just an array of stationary always-on point sources.
(Incidentally, the extra constant in the convolution kernel was chosen so that G is a Green's function for the universal wave equation: (i.e. $(\frac {\partial^2} {\partial t^2} - \nabla^2)G = \delta(x)\delta(t)$. As a result, $\ddot v - \nabla^2 v = u$.)
