Algorithm for finding initial conditions of differential equations given trajectory

Question

Let's say I'm given a system of three first-order differential equations in three variables, where all of the equations are known, and we additionally know the trajectory of two of the variables at a set of timepoints. Is there a fast algorithm to determine the starting value of the third variable, whose trajectory is unknown.

The best I've been able to come up with so far randomly initializing the starting point of the third variable, and then using a numerical solver to determine the hypothetical trajectory that the two variables with known trajectory follow, and then optimizing on the sum of the difference of squares between the hypothetical trajectory and the actual trajectories of these two variables. However, this can take a while to converge, which is why I'm looking for something faster.

Answer 1

As Christian Clason alluded to in his comment, you have reinvented the ensemble Kalman filter! The Kalman filter is an algorithm for estimating, from imperfect measurements, the true state of a system governed by a set of first-order ODEs or difference equations. For example, you might only be able to make noisy measurements of the state of the system, or you might only be able to measure some of the state variables but not all. The measurements $y$ are assumed to be a linear function of the state variables $x$, plus some noise $\eta$:

$y = Hx + \eta$

where $H$ is the measurement matrix. In your case, the measurement matrix is not of full rank.

One of the great parts about the Kalman filter is that it gives you not just an estimate of the state of the system but a probability distribution, assuming that the errors are normally distributed, even when the measurement matrix is not of full rank. While the backward ODE may be unstable, by also tracking covariance of the Kalman filter estimates of the system state, you can also put some bound on the uncertainty of your guess for the initial condition. This doesn't remedy your lack of information; it only tells you how this lack of information, i.e. only being able to observe the first two variables, results in uncertainty in the third variable.

Kaipio and Somersalo have a chapter on Kalman filtering. I also found these notes to be really helpful. If your ODE system is nonlinear, it might be beneficial to try implementing the algorithm for linear systems of ODE first. It's especially valuable to see how the results change when the governing system is stable, unstable, or oscillatory.

Answer 2

Suppose your differential equations are: $$ \dot{x} = f_x(x,y,z)\\ \dot{y} = f_y(x,y,z)\\ \dot{z} = f_z(x,y,z), $$ and you know $x$ and $y$ at (ordered) times $t_1, t_2, …$.

If $t_{i+1}-t_i$ is sufficiently small for some $i$ (preferably all of them), you have two ways to estimate the derivative of the first variable at $:

$$ \begin{alignat*}{1} \dot{x}\left( \frac{t_i+t_{i+1}}{2}\right ) &\approx \frac{x(t_{i+1})-x(t_i)}{t_{i+1}-t_i} \\ \dot{x}\left( \frac{t_i+t_{i+1}}{2}\right ) &\approx f_x \left( \frac{x(t_i)+x(t_{i+1})}{2}, \frac{y(t_i)+y(t_{i+1})}{2}, z\left( \frac{t_i+t_{i+1}}{2}\right ) \right). \end{alignat*} $$

If you equate the right-hand sides, the only unknown in this equation is $z\left( \frac{t_{i+1}+t_i}{2}\right )$, so all you have to do is to solve it for this term. You can then apply the same approach to $y$ to get a second estimate and refine your result.

(While I used the middle of the interval, you can apply the same approach to the beginning or end as well. This has less accuracy but may be more comfortable and even faster, as it saves one step in the following.)

Once you have this estimate for $z\left( \frac{t_{i+1}+t_i}{2}\right )$, you can integrate backwards and forwards to the known time points and refine your result. The best strategy here depends on the accuracy of your data and how your times $t_i$ are distributed. If $t_1$ is not your initial time, first obtain an accurate estimate of $z(t_1)$ and then integrate backwards to your initial time. This may not be much different from what you have been doing, but you should have a better initial guess.

Be aware of the possibility that your backwards differential equation may be unstable (e.g., for strongly dissipative systems), in which case you may want to avoid backwards integrations as far as possible in the above (or only perform them for very short time steps), in which you are essentially back at your method. Also, if $t_1$ is much later than your initial time (with respect to instability time scales), you have the more fundamental problem that the solution to your problem is overly sensitive on measurement and numerical errors – every solution to your problem would also pose a method to backwards-integrate unstable differential equations.